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detonating a bare charge in space. It is suggested that the loading impulse may be approximated by 
the total momentum of that portion of the fluid which impacts at the target. Assuming impulsive 
dynamic response, and assuming that the ensuing damage is proportional to the kinetic energy 
imparted to the structure by the blast, we get a particularly simple law : Damage ~ ( W is 

charge mass, R is range). This model is an idealization of a solar panel (or antenna) extended in a 
paddle-like fashion from a relatively rigid and massive core structure. It is also shown that this law 
implies that no advantage can be realized by re-arranging the mass of a single bare charge in a cluster 
configuration of smaller sub-charges, which would be dispersed and detonated via an idealized 
"isotropic" scheme. 
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1. INTRODUCTION 



The advent of space-based weapon systems in our times has raised the prospects of future 
"Star Wars" conflicts, rendering the potential use of explosive devices against space targets a present 
day engineering reality. The warhead of choice in space seems to be of the fragmentation type, for 
obvious reasons. The effectiveness of fragments is unhampered by the space environment (lack of air 
may even be helpful). By contrast, bare charges in space are considerably less efficient than in air. 
One may wonder why this is so since in air, as in space, the same amount of chemical energy is 
released through the detonation process. The explanation is that the difference is in the much larger 
mass involved in the air blast, relative to the bare charge mass. 

For a more comprehensive explanation, we take a close look at the process by which an explosive- 
driven air blast w’ave is generated. The explosive products effectively constitute a rapidly expanding 
spherical piston (typical initial speed around 6 km/sec), which drives an intense shock wave into the 
surrounding air. At a typical range of lOOR^ (and with air density equal to about 1/1000 of charge 
density), the mass of air entrained by the shock is about 1000 times the charge mass. Thus, the 
highly concentrated initial explosive energy, has spread over a much larger mass than that of the 
charge, via the mechanism of wave propagation in compressible media, resulting in an increased 
momentum. For a comprehensive treatment of blast waves in air the reader is referred to Bakerfl]. 

It is also worthwhile noting that explosive products in space typically attain hypersonic speed prior 
to impacting at the target. The flow velocity in an air blast is typically subsonic or somewhat 
supersonic. It is thus expected that the actual gasdynamic interaction between the blast flow and a 
stationary target, will be fundamentally different in these two cases. 

We contend that blast effects in space may still be of practical interest for reasons such as the 
following : 

(i) Notwithstanding the poor efficiency of a bare charge, its use should not be ruled out altogether. 
Fragments would contribute to existing - and potentially hazardous - population of space 
debris, underlining the obvious fact that there is no absolutely safe standoff distance from an 
isotropic fragmentation warhead. A clean bare charge may thus be a reasonable alternative. 

(ii) Even a fragmentation warhead has some residual blast capacity, which has to be considered 
either as a factor in enhancing target damage, or as a threat to be reckoned with in determining 
a safe standoff distance. 
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The key idea of the present model is a combination of the assumption that target dynamic 
response is related primarily to total blast impulse, and the physically plausible notion that this 
impulse is equal to the total momentum of that portion of the expanding explosive products which 
impacts at the target. The sense in which this simple notion constitutes an approximation to a proper 
gasdynamic analysis of the interaction between the fluid and the target, is clarified in Ch. 2. In that 
chapter we also present an illuminating comparison between impulsive blast loading in air and in 
space. 

In order to demonstrate the ChargeMass-Range-Damage relationship implied by our impact blast 
approximation, we chose a simple target model ; A cantilever beam with a rigid-perfectly plastic 
stress-strain relationship. It represents an extended structural element such as a solar panel or an 
antenna. We make use of studies conducted by iMentel [2] and by Bodner and Symonds [3], which 
showed that by and large, the effect of accelerating the beam impulsively was to cause a rotation 
about a plastic hinge at the point of support. The final angle of rotation is generally proportional to 
the initial kinetic energy, so that equating damage with that angle, results in damage being 
proportional to the square of the impulse imparted to the target by the blast loading. A presentation 
of this dynamic response model, including a sample case, is given in Ch. 3. 

Our ChargeMass-Range-Damage relationship may imply some far-reaching conclusions when 
applied to the analysis of a more general configuration than the single-charge/ single-target case. In 
Ch. 4 we present a simple analysis of a sub- munition configuration of N bare charges, concluding 
that it seems to have no advantage in efficiency, relative to a single charge of equal mass. Sections 5 
and 6 contain conclusions and references, correspondingly. 

We conclude the introduction by listing the main assumptions made in the present study : 

(a) Blast loading and target response are uncoupled. This is true since typically the target mass is 
much larger than the mass of that portion of the explosive products which impacts on it. 

(b) Dynamic target response is independent of specific loading time history. It depends solely on 
total (time-integrated) impulse. 

(c) The target is a panel extended as a relatively supple cantilever. It is supported by a relatively 
rigid and massive core structure. 

(d) The charge is a sphere detonated at its center. The expansion is spherically symmetric. 

(e) Target surface is normal to local flow vector. 

(f) Target orbital velocity relative to the center of the charge is negligible, compared with the 
velocity of the expanding products. 
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2. IMPACT BLAST LOADING 



Consider the expanding explosive products impacting at a target as shown in Fig. 2-1. By 
regarding the fluid as an ensemble of non-interacting particles moving at velocity U(R,t) , and by 
assuming a no-rebound normal impact at the surface, the pressure time history is given by : 

P^(t) = p(R,t)[U(R,t)]2 (2-1) 



How is this simple impact mechanism related to the actual gasdynamic interaction between the 
expanding explosive products and the target? When a target is located at a range of at least several 
charge radii, two features in the free stream of the oncoming fluid are significant : The flow is highly 
hypersonic (Mach number 20 or higher), and the static pressure is very small, which means that 
P + pU^ « pu2 . These facts were born out by a numerical computation which we performed for a 
typical high explosive characterized by the following parameters : 

Pq = 1800 (kg m’^) 

”^CJ ^ 

( 2 - 2 ) 

Dcj = 8 (m ms'*) 

Qo = Dcj-/(2(y^/ - 1)] = 4 (MJ kg'*) 

Where Qq was determined by assuming that the detonation corresponded to the CJ point on the 
explosive Hugoniot curve, and that the detonation products were an ideal gas with a specific-heat 
ratio . The spherically expanding flow was computed by integrating the Euler equations for 
isentropic flow via a high-resolution conservative finite-difference scheme [4-6]. The initial conditions 
were the self-similar flow field of a just-detonated spherical charge given by Taylor [T]. The code 
GRP with which the computation was performed is described and listed in Appendix A. 

Consider the flow at a stationary target, which begins at the moment of arrival of the expanding 
explosive products (Fig. 2-2). A qualitative description of the ensuing flow pattern is made by 
observing its evolution in time. Immediately following the initial (normal) impact, the fluid is stopped 
at the target by a backward-propagating shock wave reflected from the surface. Since the target is of 
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finite extent, the fluid between the shock and the surface is accelerated laterally, and streamlines that 
tend to curve around the target are being formed. If the oncoming flow were stationary, the flow 
field would evolve toward the familiar configuration of a detached bow-shock positioned at a 
relatively narrow standoff distance from the surface. 

Let us find the post-shock pressure in these two limiting phases. In the initial phase, the fluid is 
stopped at the target by a reflected shock (Fig. 2-3a), and in the pseudo-stationary phase (Fig. 2-3b), 
the shock is stationary. In either case we find the post-shock pressure to be given by a pressure- 
recovery expression of the form : 

P 2 = apU^ (2-3) 

Where a is a constant related to the appropriate y (assuming the expanded explosive products 
are an ideal gas). The governing equations in the reflected shock case are : 

p(U + S) = PjS 

P(U + S)2=P2 (2-4) 

P(Y + l)/(y “ 1) = P 2 (strong shock) 

Where the unknowns are p, , P, , S . 

The equations for the stationary shock case are : 
pU = p,U, 

pU2 = P 2 + P 2 U,- (2-5) 

P(Y 1)/(Y ~ 1) = P 2 (strong shock) 

Where the unknowns are p, , U 2 , P 2 • Thus, solving for a in the two cases represented by 
equations (2-4) and (2-5), we get ; 
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(2-6) 



Reflected shock a = {(y + l)/2]^ 

Stationary shock a = 2/(y + 1) 

In either case, since the gas is not dense, the efTective range of y is somewhere between 1.0 and 
1.4 , so that setting a = 1 is an approximation commensurate with the overall crudeness of the 
present impact blast model. Since the flow in the layer between the shock and the target is low 
subsonic (at least it is so away from target edges), the post-shock pressure is a reasonable substitute 
for the surface pressure. Also, a = 1 is an appropriate approximation where the flow is so 
rarefied that it is collisionless. In this limit, a = 1 corresponds to full thermal accommodation of 
re-emitted molecules from a presumably cold surface. 



The foregoing analysis constitutes a justification of the impact approximation to the surface 
pressure (2-1). Now we turn to the task of evaluating the impulse which is defined as the time- 
integrated surface pressure. Using the impact approximation (2-1), the impulse is given by : 

00 CO 



I(R) = 




■ 



= p(R,t)[U(R,t)l^dt 



(2-7) 



o 



0 



Let us introduce a Lagrange mass coordinate m which enables a transformation from the Euler 
system (R.t) to the Lagrange system (m,t). The differential relation associated with this 
transformation at constant R is : 



dm = 4JtR‘p(R,t)U(R,t)dt 



( 2 - 8 ) 



Since it is assumed that the fluid is not accelerated at any (R,t) in the range of interest for blast 
loading, the velocity U(R,t) can be regarded as function solely of the mass coordinate , so 
that U(R,t) = U(m) . Using (2-8) we are then able to cast the impact blast expression (2-7) in the 
following simple and physically appealing form ; 

I(R) = Z/47tR“ 

W (2-9) 

Z = 1 U(ni)dm 
0 

The upper limit VV in (2-9), which is consistent with the upper limit ^ in (2-7), implies that the 
total impulse is somewhat overestimated, since it contains contributions from the innermost layers of 
the explosive products that will arrive at the target as t -♦ oo . 
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The total momentum Z is thus a constant which can be evaluated for any specific explosive 
charge by numerical integration. We performed this computation with the code GRP described in 
Appendix A. In doing so for the typical explosive (2-2), we found out that the impulse (2-9) was a 
reasonable approximation at ranges as low as R = 3R^ . Furthermore, it was found that Z could 
be approximated by the maximum attainable momentum for the given charge mass and 
energy W(2 Qq)*^^ , to within about 6% . Apparently, the total momentum is not overly sensitive 
to the exact velocity distribution function U(m) , so that assuming a value of Z appropriate to the 
uniform distribution U(m) = (2 Qq)*^^ is a reasonable approximation. Thus we finally arrive at the 
following closed-form approximation for the blast impulse : 



I(R) = KW(2Q(j)'/2/4jtR2 
K = 1 



( 2 - 10 ) 



Where the coefficient K is retained in order to suggest that its value be determined more accurately 
from detailed experimental or computational data, in the event that such data become available. At 
present our best estimate is K = 1. 

There is one comparison, however, which can readily be made with available data. We refer to 
impulsive blast loading in air, such as given by Baker (Ref 1, Fig. 6.3 in the supplement). The 
comparison is conveniently made with a non-dimensional form of (2-10), which is rewritten as : 

i = I(R) [AnR^-miQ^yr-] = (R/R^)2 (2-11) 

The air blast data has to be converted to the same normalization scheme as in Eq. (2-11), before 
the comparison can be made. Considering the definition of i in (2-11) above, and the definition of 
scaled range and air blast impulse (Table 6.2 of Ref 1), this conversion is done by multiplying the 
scaled air impulse and range by the following coefficients (sea-level air is assumed) : 

Impulse Multiplier p = 3(2Y)'‘^-(4rt/3)'/^ (PJp^Q^)*/^ (p^/p(j)‘/“ = .01204 

Range Multiplier 6 = (4n/3)’'^^ = 67.06 (2-12) 

p^ = 1.3 (kg m'^) P^ = 0.1 (MPa) Y = 1-4 
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The air blast conversion was done by a small code which is given in Appendix B. The air and 
space blast impulses are shown in Fig. 2-4. We note that at ranges larger than about 10 charge radii, 
the air blast impulse is higher than the space impulse, and the gap widens as the range increases. 
This observation is consistent with the qualitative explanation given in the introduction, which 
attributed this effect to the increase in the entrained air mass at higher range. At ranges lower than 
10 charge radii, the air mass is relatively insignificant, so that one may expect the blast impulses in air 
and in space to be comparable. Indeed, the inverse-square variation of impulse with range is 
apparent for the air blast at low range. In absolute values, however, the low-range space impulse is 
higher by a factor of about 1.7. This might be interpreted as indicating that choosing K = 1/1.7 
would be the appropriate "calibration". However, we do not propose to do so, since we are not able 
to trace the various factors affecting the low-range impulse as given by Baker [1]; they may somehow 
depend on the presence of air, as well as on other parameters such as target size and equation of state 
of the explosion products. 



/ 



/ 



/ 




Figure 2-1. Impact Blast Loading 
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Figure 2-2. Shock Reflection at Impact Phase 
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(a) Initially Reflected Shock (Impact) 



(b ) Stationary Shock 



Figure 2-3. Limiting Cases of Shock Reflection 
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Figure 2-4. Impulse of Normally Reflected Blast Wave at Sea-Level and in Space 
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3. TARGET DYNAMIC RESPONSE 



For the sake of constructing representative ChargeMass- Range-Damage relations from our impact 
approximation to the blast impulse (2-10), we suggest a simple idealized structure as target model. It 
is a cantilever beam made of a metal characterized by a rigid- perfectly plastic stress-strain relation. 

This model is supposed to represent an extended spacecraft component such as a solar panel or an 
antenna. The core structure is assumed to be much more massive and rigid than the extended 
structural element, so that the cantilever can be idealized as being rigidly supported. The sole 
dynamic and structural parameters are hence those of the cantilever. 

For this purpose we make use of an experimental and theoretical investigation of uniform 
cantilever beams subjected to impulsive loading that was conducted by Mentel [2]. Aluminum alloy 
beams were held in a massive support that was gliding along a rail at speed V , until it was abruptly 
stopped by a very massive anvil. After the system came to rest, the beams were observed to have 
rotated through an angle 0 about the point of support, with little deformation elsewhere (Fig. 3-1). 

The theoretical model suggested by Mentel [21 for predicting 0(V) , can be described as 
comprising two stages. Immediately following the impact, the beam commences rotating rigidly 
about the support point, with an angular momentum equal to the pre-collision moment of 
momentum about that point. This application of the principle of conservation of moment of 
momentum entails an abrupt re-distribution of velocity in the beam, with velocity being proportional 
to distance from support, and the tip moving at 1.5 V . The angle 0 is subsequently determined 
from the requirement that the rotational kinetic energy be dissipated as plastic hinge work Mp0 . 
The resulting 0(V) expression is ; 



We now make one more step in formulating the model, in that we postulate that the angle 0 is 



0 = (3/8)pLv2/M 



P 



(3-1) 



a measure of damage . Using the following expressions for Mp , ji and V : 



Mp = (l/4)Ylr 




(3-2) 



V = I(R)/p 
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We get from (2-10) and (3-1) the following ChargeMass- Range-Damage (W-R-6) relationship : 



R = cW*/2 

C = 1(3/ 167120) (LQ^,/PpYh^)l>/4 



(3-3) 



We note that the effective range for a specified target and "damage level" 0 , is proportional to 
the square root of the charge mass W . 



Using the data for the typical explosive (2-2), and the following data for a specific aluminum beam, 
we get for this sample case : 



h = 0.002 (m) 



L = 1.0 (m) 

Pp = 2700 (kg m-2) (3-4) 

Y = 300 (MPa) 

C = 1.85 (m kg-‘/2) 

The ChargeMass-Range-Damage relationship corresponding to this sample case is depicted in 
Fig. 3-2. 
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4. CLUSTER CONFIGURATION 



In a cluster configuration, the gain in damage is presumably a result of a favorable design tradeoff 
between reduced charge mass and reduced range. Can such a gain be achieved for a space system, 
assuming the ChargeMass-Range-Damage law (3-3) to hold? It can be shown that by adopting 
some simple strategy of sub-munition dispersion and initiation, equation (3-3) implies no gain in 
target damage. 

Let us assume for the sake of a reasonably simple analysis, that dispersion and initiation of sub- 
charges would take place according to the following scheme ; 



(a) The N sub-charges appear to fan out from a common virtual center, moving at equal speeds. 
At subsequent times, their centers are uniformly distributed over an expanding spherical 
envelop. 

(b) The target moves at a constant velocity relative to the virtual center. Its point of closest 
approach to that center is at range R . 

(c) The timing for dispersion is chosen so that the target intersects (tangentially) with the spherical 

envelop at the point of closest approach (Fig. 4-1). This is also the point at which the blast 

from a single-charge configuration detonated at the virtual center, would have impacted at the 
target. 

(d) All sub-munitions are detonated at this "moment of closest approach". 

(e) It is assumed that each spherical cap of area 47tR~, N will contain one, and only one. sub- 

charge. The probability of the charge location on that cap is assumed to be uniformly 
distributed. The expected location on the cap is hence that latitude line (p which divides the 
cap into two parts of equal area (Fig. 4-2). 

(f) It is assumed that the target is subjected to the blast of a single sub-charge, which is located on 
the mid-area latitude (p of the spherical cap that surrounds the target (Fig. 4-2). 

Since the area of the spherical cap subtended by (p is 4;tR”/(2N) , the angle (p is given by : 

sin(<p/2) = (2N)-'/2 (4-1) 

We seek a comparison between the deflection 0 for a single charge (W,R) , and the deflection 

0^. in the sub-munition case ( W^. = W/N , R;>,j = 2Rsin((p/2) ). From the ChargeMass-Range- 

Damage law (3-3), using also Eq. (4-1), we get ; 
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(0^70) = (VV^/W)2 (R/R^f = 1/4 



(4-2) 



Consequently, there is no potential gain in a tradeoff between charge mass and range, for a cluster 
configuration with the aforementioned dispersion scheme. The factor 1/4 , along with the mass 
overhead inherent in constructing a multi-charge configuration, indicate that in causing blast damage, 
a single charge is more effective than an equal-mass isotropically dispersed cluster. 
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5. DISCUSSION AND CONCLUSIONS 



Our analysis pertains to a bare explosive charge initiated at a point of closest approach to the 
target. We have shown that the loading impulse on a planform target is given by the impact 
approximation (2-7), which states that the impulse is proportional to the charge mass and inversely 
proportional to the range squared. The impulse in space has been compared with impulse in air at 
sea-level. It was found that the two are quite comparable at close range (10 charge radii or less), 
exhibiting identical variation with range. At far ranges, the impulse in air is the higher one. This is 
consistent with the notion that spreading the explosive energy over larger air mass results in larger 
momentum (and hence reflected impulse). We then proceeded to develop the ChargeMass-Range- 
Damage law (3-3) for an impulse-responsive target, which states that blast damage is proportional to 
the square of the charge mass and inversely proportional to the fourth power of the range. These 
results were obtained by introducing extensive simplifications in the analysis of gasdynamic 
interaction, and in the analysis of dynamic target response. We have further shown that this damage 
law also implies that no gain can be achieved by an idealized cluster configuration of bare sub- 
charges, relative to a single charge of equal total mass. 

It is worthwhile noting that all assumptions introduced in the course of formulating the impact 
blast approximation and the structural dynamic response to impulsive loading, imply that target 
damage is overestimated. The only exception is the approximation in setting a = 1 , which can be 
readily rectified by assigning to a the reflected shock value given in (2-6). Furthermore, we assumed 
that the pressure at the midpoint of the target, is the pressure everv'where on the target. Due to flow 
around the edges, the average pressure is lower than the midpoint pressure. Also, targets are not 
everv-where normal to the flow (and charge/ target attitude is not a design parameter). Oblique impact 
obviously entails reduced target loading. In the area of structural dynamic response, a time- 
distributed loading function generally delivers less kinetic energy to the structure than an impulsive 
loading of equal total impulse, resulting in reduced deformation (damage). Thus, while the present 
model may be regarded as an over estimate when applied to a sure-fail analysis, it is particularly 
suitable in determining a sure-safe range. 
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APPENDIX A. The GRP Code 



The purpose of this Appendix is to provide a concise description of the GRP code, and a listing of 
its CHARGE version. It is intended for users that have had prior experience in implementing 
schemes for solving the Euler equation of compressible flow. The theoretical background of GRP 
schemes constitutes the principles on which the code is founded. Some familiarity (at least) with this 
background, as given in References 4, 5 and 6 is indispensable to any implementation of GRP 
schemes. Reference 4 is recommended as an introduction. The planar GRP scheme is fully described 
in Reference 5, and the duct-flow GRP scheme on which the present CHARGE version is based is 
given in Reference 6. (In CHARGE version the flow is spherical and the "duct" area is set to 
X(I)**2 , but the code can handle any area variation - see subroutines CROSS and RATIO below). 

In GRP schemes, second-order accuracy is achieved by considering a piecewise linear interpolation 
of the flow in each cell (Fig. A-1), from which second-order accurate fluxes at each cell interface are 
evaluated through an analysis of a local Generalized Riemann Problem (GRP). Briefly stated, the 
GRP goes one step further than the Riemann Problem (RP), in that it seeks (analytically) the first 
time-derivative of the flow that evolves as the "diaphragm" is removed from the cell interface, at the 
origin of the centered (X,T) wave paths of the RP solution. The major computational subroutines are 
CYCEUL where the integration of conservation laws is performed. RIEMAN where the local 
Riemann Problems are solved by Newton-Raphson iterations. MAGA where the closed-form 
expressions derived from the GRP analysis [6] are used to compute flow time-derivatives along the 
contact surface, FLL'XE where all the previously computed information is used to extrapolate the 
fluxes to mid-time-step (T + DT/2) which constitutes a second-order accurate flux. 

The plan of this Appendix is as follows. .Array variables, including those which carrv’ conserved 
variables (mass, momentum and energy), are described in section A.l. This is followed by 
descriptions of general parameters (.A. 2), labeled COMMON variables (A. 3) and all subroutines (A. 4). 
We conclude by giving the CHARGE version listing (A. 5), which should be consulted whenever a 
reading of this code description is attempted. 

NOTE : The present CHARGE version was implemented in a GRP code version that had been 
converted to treat detonation waves as chemically reactive compressible flow. However, the 
detonation scheme is effectively neutralized by setting QDET = 0 (in NETUNM). All variables 
pertaining to detonation, such as arrays Z(I), DZ(I), FIMZ(I), ZMDOT(I) and labeled CO.MMON 
variables containing Z in their names, should be ignored. 
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A.l Array Variables 



The code GRP is organized so that all major subroutines are called with standard list of array 
variables which represent the integration scheme (i.e. the conservation laws), local Riemann Problem 
solutions and second-order accurate fluxes. Virtually all array variables are initially defined in 
BEGIN (initial conditions), and are subsequently updated at each time step in CYCEUL. The 
following list explains the meaning of these variables. Some terms used in the list are defined below. 



X(I) 

U(I) 

P(I) 

RO(I) 

E(I) 



DU(I) 

DP(I) 

DRO(I) 

DG(I) 

DXSl(I) 

MlN(l) 

US(I) 



PS(1) 



LTDOT(l) 

PIDOT(I) 

FlMZ(l) 

ZMDOT(I) 



grid point coordinate, 
velocity in cell I. 

pressure in cell I (computed from equation of state). 

density in cell I. This variable is time-integrated according to the law of 
conservation of mass. (Computed in CYCEUL). 

total energy per unit volume (including kinetic energy) in cell I. This variable is 

time-integrated according to the law of conservation of (total) energy. (Computed 

in CYCEUL). 

velocity difference in cell I. 

pressure difference in cell 1. 

density difference in cell I. 

Lagrange sound velocity difference in cell I. 

the Lagrange coordinate increment defined as RO(l)*(X(I+ l)-X(I)), for cell I. 
inactive in present version. 

velocity at the contact surface obtained after the resolution of the local 

discontinuity at X(l) (Riemann Problem solution). It is denoted as U in 

References 4-6. 

pressure at the contact surface obtained after the resolution of the local 

discontinuity at X(I) (Riemann Problem solution). It is denoted as P in 

References 4-6. 

time derivative of US(I) along the contact surface. (This derivative is the result of 

the GRP analysis. It is computed in .VIAGA. See Ref 5 and 6). 

time derivative of PS(I) along the contact surface. (This derivative is the result of 

the GRP analysis. It is computed in .VIAGA. See Ref 5 and 6). 

inactive in present version. 

inactive in present version. 
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TENA(I) 

FIRO(I) 

FIM(I) 

FIE(I) 

GIP(I) 

VOL(I) 

Z(I) 

DZ(I) 



momentum per unit volume RO(I)*U(I) in cell I. This variable is time-integrated 
according to the law of conservation of momentum. (Computed in CYCEUL). 
mass flux at point X(I) (second-order accurate), 
momentum flux at point X(I) (second-order accurate), 
energy flux at point X(I) (second-order accurate). 

the pressure term in the momentum flux. It corresponds to G(U) in References 4 
and 6. 

volume of cell I. 
inactive in present version, 
inactive in present version. 



Glossary of terms used in the array variables list : 

Cell I - the cell between grid points X(I) and X(I + 1). All cell variables are averages per that 
interval. 

Difference in cell I - the difference between values of variable at cell boundaries X(l + 1) and X(I). 
Those values are obtained from "monotonized" piecewise linear distribution of each variable 
in each cell. (Fig. A-1). 

Second-order accurate flux - the flux time-derivative at point X(I) is computed from the time- 
derivatives of pressure and velocity along contact surface PIDOT(I) and LTDOT{I ) (in 
FLUXE). Then the the flux is extrapolated to the centered time point (T + DT 2). using 
those derivatives. This centered value is the second-order flux for integrating the 
conservation laws between T and T + DT. 
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A.2 Major Parameters 

A list of major parameters indicating their meaning and the routine in which they are defined, is 
given below. Those parameters defined in KETUNM are the run input. There is no reading of an 
input file in this version of GRP code (and the only output is the printed output). 



L 

LL 

T 

DT 

TMAX 

TMUD 

DTMUD 

XCYC 

COLELA 

KEYMON 

NCYCPR 

STAB 

DTBA 

DTKOD 

KDT 



number of grid points + 1 (main program) 

L - 1 (MAIN PROGRAM) 
time (MAINO) 
time step (MAINO) 

maximum time (when T.GE.TMAX the run is terminated) (NETUNM) 
time for which next printing will take place (NETUNM) 
printing time step (NETUNM) 

serial number of time step (integration cycles) (MAINO) 

switch to evaluate cell differences by Colella's method when COLELA.NE.O 
(NETUNM) 

key for monotonization scheme (just one is presently provided when COLELA. EQ.O) 
(NETUN.M) 

frequency of line printing at each cycle (time step) (NETUNM) 

CFL stability coefficient. Must be smaller than 1. (NETUN.VI) 
next time step computed from stability criterion (CYCEUL) 
former time step (MAINO) 

index of cell where DTBA was determined (CYCEUL) 
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A.3 Labeled COMMON variables 



Labeled COMMONS are used primarily to transmit data to and from routines that perform the 
major computational steps of the GRP scheme, i.e, RIEMAN, MAGA and FLUXE; these routines 
are called from CYCEUL. When the value of any of those variables is needed for later use, whether 
for updating conservation variables (RO, TEN'A, E), or for printing, it is stored in the appropriate 
array. All labeled COMMON variables are grouped under labels that indicate their role, and their 
names are also mnemonic. Generally, sufTi.x L means Left and suffix R means Right. It may indicate 
sides either with respect to a cell interface X(I), or with respect to the contact surface which separates 
the Right- and Left- propagating weaves in a solution to the local Riemann Problem. We indicate by 
INPUT variables that are computed prior to calling the subroutine, and by OUTPUT variables whose 
value w^as computed within the subroutine and constitutes the result of calling that subroutine. 



COMMON /STEPO/ Parameters related to the local Riemann Problem. This is the first step in 
the GRP scheme. 

UL, PL, ROL, CL, GL, SL - velocity, pressure, density, sound speed, Lagrange sound speed and 
entropy, attributed to Left side of cell interface at point X(l). (INPUT) 

USTAR, PSTAR - velocity and pressure at the contact surface obtained when the local 
discontinuity is resolved (i.e.. the solution to the local Riemann Problem). The 
omission of L or R sutfix indicates that P and U are continuous across the contact 
surface. (OUTPUT) 

RSTARL, CSTARL, GSTARL - density, sound speed and Lagrange sound speed on the Left side of 
the contact surface. (OUTPUT) 

WL - Lagrange velocity of propagation of the Left-moving shock, relative to the fluid. (OUTPUT) 
UW(6) - velocity of propagation of each wave front (Fig. A-3), relative to the inertial system 
(X). (OUTPUT) 

HELEML - logical variable. If HELEML.EQ..TRUE. the Left-propagating wave is a shock. 

Otherwise it is a (centered) rarefaction wave. (OUTPUT) 

NFLUX - integer variable. It denotes the region in the Riemann solution wave structure, which 
contains the point X(I) for all time. Refer to Fig. A-3 for illustration. (OUTPUT) 
LAMDAL. R.ATEL, TEMPL, TEMPSL. ZL, ZSTARL - inactive. 
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COMMON /STEPl/ Parameters related to the time-derivative evaluation of the GRP scheme, 
performed in MAGA. The time-derivatives of P and U along the contact surface are 
the main result of MAGA. 

DUIDT, DPIDT - time-derivatives of velocity and pressure along contact surface. (OUTPUT) 

ASTARL - The directional derivative of U along the fan characteristic at the trailing characteristic 
of the Left rarefaction wave. It is not evaluated when the Left wave is a shock. (See 
References 4-6) (OUTPUT) 

DGIDTL, DRIDTL - time-derivatives of Lagrange sound speed and density along the left side of 
the contact surface. (OUTPUT) 

DSDAL - Lagrange spatial derivative of entropy on the left side of contact surface, prior to 
removal of the partition at X(I). 

SH, RAT - the cross-section area and the x-derivative of In(SH). They are user-defined in CROSS 
and RATIO respectively. 

DSDASL - entropy derivative used in the special "sonic" case (i.e, when NFLUX = 2 or 
NFLUX = 5). See References 5,6 for details. (OUTPUT) 

LAMDSL, DZDAL, BETACL, DZDASL - inactive. 



COMMON /GRADS/ Used to transmit flow gradients (that exist in fluid prior to removal of the 
partition at X(I)) to M.A.GA. 

DUDXIL, DPDXIL, DGDXIL, DRDXIL. DSDXIL - gradients of U. P. G. RO. S (with respect to 
Lagrange coordinate). They are computed in CYCEUL for transmission to .VIAGA. 
(INPUT) 

DZDXIL - inactive. 



COMMON /FI/ Used to return values of updated flux and cell-interface variables from 

FLUXE. 

FIHl, FIH2, FIH3 - second-order flux of mass, momentum flow (just RO*U’'”‘‘2) and energy. They 
are extrapolated to Half the time step T + DT/2. (OUTPUT) 

GIH - the value of P at T+ DT/2 
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UXN, PXN, GXN, ROXN - values of U, P, G, RO extrapolated to New time T + DT, at cell- 
interface. They are used in CYCEUL to get tentative (pre-monotonized) new cell 
differences. .(OUTPUT) 

ZXN, FIH4, ZMDOTL, ZMDOTR - inactive. 
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A.4 Description of Subroutines 



MAIN PROGRAM 

The task of this program is to allocate array space for the NMAT arrays required by the present 
version of GRP code. The length of each array is L. The allocation is done by calling MAINO. This 
standard calling sequence is maintained hereafter, thus facilitating modifications. 

MAINO 

This subroutine functions as an overall organization routine. It can be read as a kind of flow- 
chart of the entire computation. First, run set-up is done by calling once to NETUNM (data) and 
BEGIN (initial conditions). Then a loop over time steps is begun. In each cycle the integration by 
one time step is performed by calling CYCEUL, and subsequently boundary conditions are 
implemented by calling SAFAE. Whenever T.EQ.TMUD, results are printed by calling PRINT and 
TMUD is updated by adding DTMUD. 

NETUNM 

Here data are set for a particular run. User is invited to modify this routine. There is no input 
file. This routine is called just once from MAINO. Note that the detonation data section is skipped 
when QDET.EQ.O. 

BEGIN 

Initial conditions are set-up in this routine. The configuration of some nominal case is given in 
present version. (In CHARGE version it is the detonated spherical charge, using the Taylor self 
similar solution as initial conditions). User is called to modify this routine so as to generate any other 
desired initial configuration. 

TAYLOR 

The purpose of this routine, along with ancillary routines INIDAT, RUNGE and DERIV, is to 
compute the self-similar Taylor solution [7] of a detonated spherical charge, and implement it as 
initial conditions for the GRP computation of the ensuing expansion. TAYLOR is called once by 
BEGIN. 

The core of the solution is the numerical (Runge-Kutta) integration of two coupled ordinary' 
dilTerential equations. The integration variable is PSI. (The flow velocity normalized by DCJ is given 
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by U = EXP(-PSI) ). The two dependent variables are X - the normalized radial coordinate (X= 1 at 
the sphere boundar>'), and C - the normalized speed of sound. The integration is carried out by 
calling RUNGE, which in turn calls DERIV for the evaluation of derivatives. Data for the TAYLOR 
computation is set up by calling (just once) I NT DAT. 

The initial conditions needed in BEGIN are values of mass, momentum and (total) energy per cell. 
These are most accurately computed by spatially integrating the Taylor solution, resulting in lumped 
mass, momentum and energy per cell, which are then divided by the cell volume. This refinement is 
significant since gradients are high near the charge boundary (X= 1). A total mass and energy check 
for the entire sphere is performed and printed. 

INIDAT, RUNGE, DERIV 

Subroutines used only in conjunction with the Taylor initial conditions setup. See TAYLOR 
above. 

RATIO, CROSS 

User-defined routines. If A(X) is the duct cross-section area, then CROSS(X) = A(X) and 
RATIO(X)=D[ln(A(X))]/DX . 

CYCEUL 

This is the central computation routine. All major stages of the GRP scheme are performed by 
calling specific subroutines from CYCEUL. Then RO(I), TENA(I) and E(I) are updated to new time 
T + DT by solving the appropriate conservation laws in CYCEUL. 

The first loop (DO 1) performs a set of preparatory' steps as follows : 

(a) CALL RIE.VIAN - Solving the local Riemann Problem at each X(l). 

(b) CALL M.AGA - Solving the local Generalized Riemann Problem at each X(l). 

(c) CALL FLUXE - Computing second-order fluxes at X(I). 

(d) Evaluation of cell-interface finite differences DU'(l). DP(I). DRO(I) in each cell. These will be 
used at the future time step (after monotonization) for piecewise-linear interpolation of the How 

in each cell. (See definition of DUDXIL. DPDXIL just preceding the call to MAGA in this 

loop). 

Note that in present CHARGE version additional computation of PRESS. PU’LSEl,..., PULSE4 
has been added. It is just informative and does not interfere in any way with the execution of the 
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GRP scheme. The purpose of this computation is to monitor the numerical solution and to observe 
the accuracy within which the asymptotic value of the momentum integral Z (Eq. 2-9 above) is 
approached. 

In the second loop (DO 2), the integration of the three conservation laws is performed, using 
second-order fluxes that had been computed in loop 1. Flow variables such as P(I) and U(I) are 
computed in this loop from the conserved variables. The cycle computation is concluded by calling 
BDOKl for monotonization of DL’(I), DP(1) and DRO(I). 

SAFAE 

In this routine user-defined boundary conditions are implemented. Present version (CHARGE) 
contains rigid wall at the center of the sphere X(2)=0, and an "open boundary" at the outer 
computational zone limit X(L). The rigid wall condition is achieved by setting up a virtual 
antisymmetric cell next to the boundary cell, so that the solution to the local Riemann Problem will 
result in a non-moving contact surface (USTAR = 0). The open boundary' is an approximation to an 
ideally non-reflecting boundary. Here the virtual cell is I = L, and the flow in it is defined as a 
"continuation" of the flow in the adjacent last cell I = LL. 

BDOKl 

Here the tentative cell-interface differences DV(I) are monotonized according to neighboring 
average cell values V(I-l), V(I) and V(l-f-l). The basic idea is that the cell-interface slope DV(I) 
should have the same sign as the average slope V(I-t- l)-V(I-l). When V(I) is a local extremum DV(I) 
is set to zero. Also, the absolute value of DV(I) is constrained so that the jump from a cell-interface 
value to the adjacent average value V(I), will never be of opposite sign to DV(I). 

DCOLE 

When COLELA option is used (not in present CHARGE version), the pre-monotonized slopes 
are simply the centered difference (V(I4- l)-V(l-l))/2, Note that even under this option, the 
monotonization routine BDOKl is subsequently called. 

PRINT 

Printing of results. Reading this routine is self-explanatory. Note some features added for 
present CHARGE version. User is called to modify this routine to his specific needs. 
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SOF 



Run termination when an error has been detected. ISTOP is an informative index. All printing 
of relevant information should be done at the calling routine prior to calling SOF. Note that the run 
is ended in SOF by deliberately causing a system error of computing SQRT(-l). This is done in order 
to trigger printing of the sequence of calling routines by the operating system. 



RIEMAN 



Here a single Riemann Problem (RP) is solved by calling RIEMAN from CYCEUL. Referring 
to Fig. A-2, the RP is solved by finding the point of intersection (USTAR,PSTAR) of Left- 
propagating and Right-propagating shock/rarefaction adiabats in the (U,P) plane. Prior to the actual 
computation, the qualitative wave structure is determined. It is characterized by the index NCASE as 
follows : 



NCASE =1 
NCASE = 2 
NCASE =3 



NCASE =4 



Left wave is rarefaction, Right wave is shock. 
Both waves are shock. 

Left wave is shock. Right wave is rarefaction. 
Both waves are rarefaction. 



The computation of (USTAR,PSTAR) is coded separately for each case. Newton-Raphson 
iteration is employed, the first guess being the intersection of the Left and Right rarefaction branches 
{or their extrapolations), which is done in closed-form. Since in a smooth How this guess is close to 
the exact (USTAR,PSTAR), little extra CPU effort is spent on subsequent Newton-Raphson 
iterations. These are truly needed only in regions of shock wave computation. 

The computation in RIEMAN is concluded by computing UW(1) UW(5) (UW(6)= infinity). 

From these wave speeds, the flux index NFLUX that denotes the location of the X-axis on the (X.T) 
wave diagram of the RP solution (Fig. A-3), is evaluated. It is later needed in subroutine FLL'XE. 
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MAGA 



The major purpose of this routine is to compute DUIDT and DPIDT along the contact surface 
of the RP solution. Since U and P are continuous across the contact, so are their time-derivatives 
along the contact. Thus, DUIDT and DPIDT are solved from a set of two linear equations. The 
coefficients of each equation are determined by GRP analysis of the wave on one side. See 
References 4-6 (particularly Ref 6) for details. 

FLUXE 

The major task of this routine is to compute second-order fluxes. This is done in tw'o phases. 
The first phase is up to statement 9 CONTINUE, where using NFLUX the X-suffixed values of flow 
variables and their time-derivatives are defined. An X-suffix means that the variable or its time- 
derivative are related to the line X = X(I) on the (X,T) wave diagram (Fig. A-3). In the second phase, 
these variables and their time-derivatives are used to extend fluxes at X(I) to Half-time-step (hence 
the suffix H), i.e. T + DT/2. It is these fluxes which are the second-order accurate fluxes for the 
integration of the conservation laws from T to T+DT. Also, cell-interface flow variables (suffix N) 
are extended to New time level T + DT. These are later used in defining cell differences DU(I), DP(I) 
and DRO(I) in CYCEUL. 



28 



A.5 Listing of GRP Code 



C$0PTI0NS LIST i.£Rs,o<i 

IMPLICIT REAL5(8(A-H,0-Z,$) vbK5ioji_ 

C PROGRAM GRP - GENERALIZED RIEMANN PROBLEM. 

C EXPANSION OF A DETONATED SPHERICAL CHARGE IN VACUUM. 

C INITIAL CONDITIONS FROM TAYLOR'S SELF SIMILAR SOLUTION. 

COMMON B(102,26) 

1 ,ENDB 

COMMON /AB/A(50) 

EQUIVALENCE ( L , A( 1) ) , ( LL , A(2) ) , (T, A( 3) ) , ( DT, A(4) ) , (TMAX, A( 5) ) , 

1 (TMUD,A(6)), (DTMUD,A(7)),(J0B,A(8)),(NERI,A(9)), 

2 (JJJ,A(I0)), (KEYM0N,A(1D), (NCYC, A(12)) 
(C0LELA,A(13)) 

(LAGEUL,A(14)) 

(UGAL,A(15)) 

(KEYEK,A(16)) 

(NCYCPR,A(17)) 

(STAB,A(18)),(DTBA,A(19)),(DTKOD,A(20)),(KDT,A(2D) 
/M0NIT/CASEAV(4),NC14(4),NF16(6), 

1 NM0NU(4),NM0NP(4),NM0NR0(4),NM0NZ(4) 

DIMENSION NZER0(26) 



EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
COMMON 



EQUIVALENCE ( NZERO( 1 ) , NC14( 1) ) 

COMMON/PULS/PRESS(10),PULSE1(10),PULSE2(10),PULSE3(10),PULSE4(10) 



CHAOOOl 

CHA0002 

CHA0003 

CHA0004 

CHA0005 

CHA0006 

CHA0007 

CHA0008 

CHA0009 

CHAOOlO 

CHAOOll 

CHA0012 

CHA0013 

CHA0014 

CHA0015 

CHA0016 

CHA0017 

CHA0018 

CHA0019 

CHA0020 

CHA0021 

CHA0022 



20 

21 
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DO 20 N=l,26 
NZERO(N)=0 
DO 21 N=l,4 
CASEAV(N)=0. 

DO 31 N=l,10 
PRESS(N)=0 . 

PULSE1(N)=0 . 

PULSE2(N)=0. 

PULSE3(N)=0. 

PULSE4(N)=0. 

CONTINUE 

NMAT=26 

L=(LOCF(ENDB)-LOCF(B(I,l)))/NMAT 

L = 102 

LL=L-1 

NN = NMAT5«L 

DO 1 1 = 1, L 

DO 1 11=1, NMAT 

B(I,II)=0. 



CALL 


MAIN0(L,B(1 




1) 


,B(1 


, 2) 


,B(] 


L, 3) 


, BC] 


L. 4) 


,B(] 


L, 


5 


1 


B(l, 


6 


), 


B(l, 


7), 


B(l, 


. 8), 


B(l; 


' 9), 


B(l, 


.10 


) 


2 


B(l, 


11 


), 


B(1 . 


12), 


B(l, 


.13), 


B(l; 


.14), 


B(l; 


.15 


) 


3 


B(l, 


16 


), 


B( 1, 


17 ), 


B(l, 


.18), 


B(1 ; 


.19), 


B(l; 


.20 


) 


4 


B(l, 


21 


), 


B(l, 


22), 


B(l, 


.23), 


B(l; 


.24), 


B(l; 


.25 


) 


5 


B( 1, 


26 


)) 





















STOP 

EN D 



CHA0024 

CHA0025 

CHA0026 

CHA0027 

CHA0028 

CHA0029 

CHA0030 

CHA0031 

CHA0032 

CHA0033 

CHA0034 

CHA0035 

CHA0036 

CHA0037 

CHA0038 

CHA0039 

CHA0040 

CHA0041 

CHA0042 

CHA0043 

CHA0044 

CHA0045 

CHA0046 

CHA0047 

CHA0048 

CHA0049 

CHAQQ5Q 



SUBROUTINE MAINO 

1 (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN, 

2 US,PS,UIDOT,PIDOT, 

X FIMZ,ZMDOT, 

3 TENA, FIRO,FIM, FIE,GIP, VOL,Z,DZ) 

IMPLICIT REALX8(A-H,0-Z,$) 

DIMENSION X(L),U(L),P(L),RO(L),G(L),E(U,DU(L),DP(L),DRO(U, 

1 DG(L),DXSI(L),MIN(L), 

2 US(L),PS(L),UIDOT(L),PIDOT(U 

3 ,TENA(L),FIRO(L),FIM(L),FIE(U 

4 ,GIP(L),VOL(L),Z(L),DZ(L) 

5 , FIMZ(L),ZMDOT(L) 

COMMON /AB/AC50) 

EQUIVALENCE 



(LL, A(2) ), (T, A(3) ), (DT,A(4) ), (TMAX,A(5) ), 

1 (TMUD,A(6)),(DTMUD,A(7)),(JOB,A(8)),(NERI,A(9)), 

2 (JJJ,A(10)),(KEYMON,A(11)),(NCYC,A(12)) 

(LAGEUL,A(14)) 

(NCYCPR,A(17)) 

(STAB,A(18)), (DTBA,A(19)), (DTKOD,A(20)),(KDT,A(2D) 

COMMON /TOT/AMTOT, ETOT, EKTOT, EPTOT,TENTOT 

T=0. CHA0072 



EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 



CHA0051 

CHA0052 

CHA0055 

CHA0054 

CHA0055 

CHA0056 

CHA0057 

CHA0058 

CHA0059 

CHA0060 

CHA0061 

CHA0062 

CHA0065 

CHA0064 

CHA0065 

CHA0066 

CHA0067 

CHA0068 

CHA0069 

CHA0070 
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FILE: 



CHARGER FORTRAN A1 



NCYC=0 
JJJ = 0 

CALL NETUNM 
DELT=DT 
CALL BEGIN 

1 
2 
X 

3 

CALL SAFAE 

1 

2 
X 

3 

I NCYC=NCYC+1 
C TIME STEP CONTROL. 

DT=DTBA 

IFCDT.GT.l .IDOXDTKOD.AND.DTKOD.NE.O.) DT=1 . IDOXDTKOD 
IF(NCYC.EQ.2) DT=DT/10.D0 
IF (NCYC.EQ.l) DT=0. 

IFCDT.EQ.O.) GO TO 11 
NHAD=((TMUD-T)/DT-1.D-10) 

IF(NHAD.GE.IO) GO TO 11 
DT=(TMUD-T)/DFLOAT(NHAD+l) 

II CONTINUE 
T=T+DT 

IF((NCYC/NCYCPR)xnCYCPR.NE.NCYC.AND.NCYC.GT.NCYCPR) GO TO 33 
PRINT 10, NCYC,T,DT,KDT 

10 FORMATCIX, ' NCYC= , 3X, *T= • , Dll . , 3X, » DT= • , Dll . 4, 3X, •KDT=M4) 
33 CONTINUE 

DTBA=DTMUD 

KDT=0 



(L,X,U,P,R0,G,E,DU,DP,DR0,DG,DXSI,MIN, 

US,PS,UIDOT,PIDOT, 

FIMZ,ZMDOT, 

TENA,FIR0,FIM,FIE,GIP,V0L,Z,DZ) 

(L,X,U,P,R0,G,E,DU,DP,DR0,DG,DXSI,MIN, 

US,PS,UID0T,PID0T, 

FIMZ,ZMDOT, 

TENA,FIR0,FIM,FIE,GIP,V0L,Z,DZ) 



NERI=1 

IF (DABS(T-TMUD) .LT.l .D-8) NERI=0 
CALL CYCEUL 



1 

2 

X 

3 

CALL SAFAE 



(L,X,U,P,R0,G,E,DU,DP,DR0,DG,DXSI,MIN 

US,PS,UIDOT,PIDOT, 

FIMZ,ZMDOT, 

TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ) 



1 

2 

X 

3 

IF (NERI.NE 
CALL PRINT 



(L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN, 

US,PS,UIDOT,PIDOT, 

FIMZ,ZMDOT, 

TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ) 

.0) GO TO 2 



1 (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN, 

2 US,PS,UIDOT,PIDOT, 

X FIMZ,ZMDOT, 

3 TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ) 

IF (DABS(T-TMUD) .LT. 1 .D-8) TMUD=TMUD+DTMUD 
CONTINUE 

DTKOD=DT 

IF (T.LT.TMAX-1 .D-8) GO TO 1 
RETURN 
END 






SUBI^OUTINE NETUNM 
IMPLICIT REAL3«8(A-H,0-Z,$) 

COMMON /AB/A(50) 

EQUIVALENCE (L,A(D) 

EQUIVALENCE ( LL , A( 2) ) , (T, A( 3 ) ) , ( DT, A( 4 ) ) , ( TMAX, A( 5) ) , 

1 (TMUD,A(6)),(DTMUD,A(7)), (J0B,A(8)),(NERI,A(9)), 

2 (JJJ,A(10) ), (KEYMON,A(ll) ), (NCYC,A(12) ) 

EQUIVALENCE (C0LELA,A(13) ) 

EQUIVALENCE ( LAGEUL , A( 14) ) 

EQUIVALENCE (KEYEK, A( 16 ) ) 

EQUIVALENCE (NCYCPR, A( 17 ) ) 

EQUIVALENCE ( STAB , A( 18 ) ) , ( DTBA, A( 19 ) ) , ( DTKOD, A( 20 ) ) , ( KDT, A( 21) ) 
COMMON/DETO/QDET,PCJDET,RCJDET,UCJDET,DCJDET,PODET,ROODET, 

1 RATE,TEMPC 

C0MM0N/DIFFUS/U2,P2,R02,ARW 

COMMON /DRAH/GODELX,GODELY, UMIN, UMAX, PMIN, PMAX, ROMIN, ROMAX 



CHA0073 

CHA0074 

CHA0075 

CHA0076 

CHA0077 

CHA0078 

CHA0079 

CHA0080 

CHA0081 

CHA0082 

CHA0083 

CHA0084 

CHA0085 

CHA0086 

CHA0087 

CHA0088 

CHA0089 

CHA0090 

CHA0091 

CHA0092 

CHA0093 

CHA0094 

CHA0095 

CHA0096 

CHA0097 

CHA0098 

CHA0099 

CHAOlOO 

CHAOlOl 

CHA0102 

CHA0103 

CHA0104 

CHA0105 

CHA0106 

CHA0107 

CHA0108 

CHA0109 

CHAOllO 

CHAOlll 

CHA0112 

CHA0113 

CHA0114 

CHA0115 

CHA0116 

CHA0117 

CHA0118 

CHA0119 

CHA0120 

CHA0121 

CHA0122 

CHA0130 

CHA0131 

CHA0132 

CHA0133 

CHA0134 

CHA0I35 

CHA0136 

CHA0137 

CHA0138 

CHA0139 

CHA0140 

CHA0141 

CHA0142 

CHA0143 

CHA0144 

CHA0145 

CHA0146 

CHA0147 

CHA0148 

CHA0149 

CHA0150 

CHA0151 
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1 

1 

2 



1 

2 

3 

4 



,XMIN,XMAX,SMIN,SMAX, IVERSA 

COMMON /GAM/GAMA , NG, MU2, G1 , G2 , G3 , G4 , G5 , G6 , G7 , G8 , G9 , G1 0 , G1 1 
,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23 
,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,G35 
REALX8 NG,MU2 

NAMELIST /IN/LIN,GAMA,DT,TMUD,DTMUD,TMAX, 

GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX, 
SMIN,SMAX, IVERSA, KEYMON,COLELA, STAB 
,LAGEUL,KEYEK 
,QDET 



CHA0152 

CHA0153 

CHA0154 

CHA0155 

CHA0156 

CHA0157 

CHA0158 

CHA0159 

CHA0160 

CHA0161 






) 



LIN=L 

LAGEUL=2 

NCYCPR=1 

KEYEK=1 

TMUD=0. 

DTMUD=10.D0 
TMAX=100 .DO 
STAB=0.5D0 
DT=l,D-2 
KEYMON=l 
GAMA=3.D0+l.D-6 
QDET=0. 04D0 
RATE=0 . 

TEMPC=1.D50 

GODELX=16DO 

GODELY=20.D0 

IVERSA=100 

UMIN=0, 

UMAX= l.DO 
PMIN=0. 

PMAX=0.5D0 

ROMIN=0. 

R0MAX=3.D0 

SMIN=0. 

SMAX=0 .03D0 
C0LELA=0 . 

READ IN 

PRINT IN 

GG=2.D0XGAMA/(GAMA-1 .DO) 

NG=GG 

CONTINUE 

MU2=(GAMA-1 .DO)/(GAMA+l .DO) 

G1=GAMA-1 .DO 
G2=l .D0-MU2 

G3 = 2.D0/(3.D05<GAMA-1 .DO) 

G4=(GAMA+1 .DO)/2.DO 

G5 = 0 . 5D05«( 3 . D05(GAMA-1 . DO )/(GAMA+l . DO 
G6=(GAMA+1.D0)/(2.D0XGAMA) 
G7=2.D0/(GAMA-1 .DO) 

G8 = ( GAMA- 1 . DO )/ ( 2 . DOXGAMA ) 

G9 = (GAMA+1 .DO)/(4. D05(GAMA) 

G10=l. DO/GAMA 
G11=(GAMA+1.D0)/4.D0 
G12=GAMA/(GAMA-1 .DO) 
G13=0.5D0X(GAMA-3.D0)/(GAMA+1 .DO) 

G14 = 0.5D0>((3.D0XGAMA-5.D0)/(GAMA+1 .D 
G15=GAMAX( 3 . DOXGAMA-1 . DO ) 

G16 = (GAMA+1 .D0)/(2.D03((GAMA-1 .DO)) 
G17=GAMA+1 .DO 

G18=GAMA5<(GAMA+1 . DO )/( 3 . D05(GAMA-1 . DO 
G19 = (3.D05<GAMA-1.D0)/(GAMA+1 .DO) 

G20 = 2.D05((GAMA-1 .D0)/(3.D0XGAMA-1 .DO 
G21=GAMA3«( 3 . DO>(GAMA-5 . DO )/( 3 . DOXGAMA 
G0DELX=G0DELX/2. 54DO 
GODELY=GODELY/2.54DO 
CALL NAMPLT(IVERSA) 

CALL LIMIT(IOOO.DO) 

CALL PLOT(0.,0.5D0,-3) 

P0DET=0. 



CHA0163 

CHA0164 

CHA0165 

CHA0166 

CHA0167 

CHA0168 

CHA0169 

CHA0170 

CHA0171 

CHA0172 

CHA0173 

CHA0174 

CHA0175 

CHA0176 

CHA0177 

CHA0178 

CHA0179 

CHA0180 

CHA0181 

CHA0182 

CHA0183 

CHA0184 

CHA0185 

CHA0186 

CHA0187 

CHA0188 

CHA0189 

CHA0190 

CHA0191 

CHA0192 

CHA0193 

CHA0194 

CHA0195 

CHA0196 

CHA0197 

CHA0198 

CHA0199 

CHA0200 

CHA0201 

CHA0202 

CHA0203 

CHA0204 

CHA0205 

CHA0206 

CHA0207 

CHA0208 

CHA0209 

) CHA0210 

CHA0211 

CHA0212 

CHA0215 

CHA0214 

CHA0215 

XXZ CHA0216 

1.D0)5«X2 CHA0217 

CHA0218 

CHA0219 

CHA0220 

CHA0221 

CHA0222 

CHA0223 
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101 

102 

103 

104 
100 



R00DET=0. CHA0224 

PCJDET=0. CHA0225 

UCJDET=0. CHA0226 

DCJDET=0. CHA0227 

RCJDET=0. CHA0228 

IF(QDET.LE.O.) GO TO 100 CHA0229 

DETONATION DATA CHA0230 

QDET=0.04D0 CHA0231 

P0DET=0. CHA0232 

ROODET=1.8DO CHA0233 

PCJDET=P0DET-(GAMA-1 . DO)X(-QDET)XROODET+ CHA0234 

1 DSQRT(( (GAMA-1. DO)XQDETXROODET)XX2-2.DOXMU2*GAMAX CHA0235 

2 (-QDET)>(PODETXROODET) CHA0236 

RCJDET=R00DETX((GAMA+1 .DO)XPCJDET-PODET)/(GAMAXPCJDET) CHA0237 

CCJ = DSQRT ( GAMA5CPCJ DET/RC J DET ) CHA0238 

DCJDET=CCJXRCJDET/ROODET CHA0239 

UCJDET=DCJDET-CCJ CHA0240 

PRINT 101 CHA0241 

FORMATCIHI,/, IX, 'DETONATION DATA'/) CHA0242 

PRINT 102, QDET, GAMA, TEMPO, RATE CHA0243 

F0RMAT(/1X, ' QDET,GAMA, TEMP, RATE= ' , 4D18 .8) CHA0244 

PRINT 103, ROODET,PODET CHA0245 

F0RMAT(/1X, 'UNBURNED STATE ROODET, PODET= ' , 2D18 . 8) CHA0246 

PRINT 104, DCJDET,PCJDET,RCJDET,UCJDET CHA0247 

FORMATC/IX, 'CJ POINT DCJDET, PCJDET, RCJDET, UCJDET= ' , 4D18 . 8 ) CHA0248 

CONTINUE CHA0249 

RETURN CHA0250 

END CHA0251 






SUBROUTINE BEGIN 

1 (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN, 

2 US,PS,UIDOT,PIDOT, 

X FIMZ,ZMDOT, 

3 TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ) 

IMPLICIT REALX8(A-H,0-Z,$) 

DIMENSION X( L ) , U( L ) , P( L) , RO ( L) , G( L ) , E( L ) , DU( L ) , DP( L ) , DRO ( L) , 

1 DG(L),DXSI(L),MIN(L), 

2 US(L),PS(L),UIDOT(L),PIDOT(L) 

3 ,TENA(L),FIRO(L) ,FIM(L),FIE(L) 

4 ,GIP(L),VOL(L),Z(L),DZ(L) 

5 , FIMZ(L),ZMDOT(L) 

/AB/A(50) 

/GAM/GAMA, NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11 
,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23 
,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,G35 
NG,MU2 

(LL,A(2)) 

(LAGEUL,A(14)) 

(UGAL,A(15)) 

(STAB,A(18)), (DTBA,A(19)) , (DTKOD,A(20)),(KDT,A(2D) 
COMMON/DETO/QDET, PCJDET, RCJDET, UCJ DET, DCJDET, PODET, ROODET, 

1 RATE, TEMPO 

COMMON /DRAW/GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX 
1 ,XMIN,XMAX,SMIN,SMAX,IVERSA 

COMMON/GIT/ROLIM,ELIM,XGIT(200),ROGIT(200),ROUGIT(200),EGIT(200) 
COMMON /GITN/NPO 
LOGICAL CSOF 



COMMON 

COMMON 



REALX8 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 



DTBA=0 . 

DTK0D=0. 

KDT = 0 
P0=1 . D-9 
RH00=1 .D-7 
UO=UCJDET 
UGAL=0 . 

X0 = 0 . 

X1=50.D0 

XCHARG=10.DO 

XMIN=XO 

XMAy=yi 

DX=(X1-X0)/(L-2.D0) 
DO 1 1=2, L 
X(I)=X0+(I-2.D0)xDX 



CHA0281 

CHA0282 

CHA0283 

CHA0284 

CHA0285 

CHA0286 

CHA0287 

CHA0288 

CHA0289 

CHA0290 

CHA0291 

CHA0292 

CHA0293 

CHA0294 

CHA0295 
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;ep fortran ai 



CONTINUE 

X(L)=X1 

UO=UOXCROSS(XCHARG) 

DO 2 1=2, LL 

IF(I.GT.2) U(I)=UO/CROSS(X(D) 

P(I)=PO 

RO(I)=RHOO 

Z(I)=0. 

GO TO (31,32), LAGEUL 
;i CONTINUE 

E(I)=P(I)/((GAMA-1.D0)xRO(I))+0.5D0XU(I)xx2+Z(I)xQDET 
GO TO 30 
12 CONTINUE 

E(I)=P(I)/(GAMA-1 .D0)+0 .5DOXRO(I)XU(I)XX2+Z(I)XRO(I)XQDET 
10 CONTINUE 

G(I)=DSQRT(GAMAXp(I)XRO(D) 

; CONTINUE 

DO 3 1=2, LL 
TENA(I) = RO(I)3(U(I) 

VOL(I)=(X(I+l)-X(I))x(X(I+l)xx2+X(I+l)xX(I)+X(I)xx2)/3.DO 
1 CONTINUE 

INSERT DETONATED CHARGE FLOW FIELD FROM TAYLOR'S SOLUTION. 

CALL TAYLOR(GAMA) 

RONORM=RCJDET 

RUN0RM=RCJDETXUCJDET 

ENORM=RCJDET3«DCJDETXX2 

XLIM=XGIT(NPO) 

NGIT=NP0-1 

XG1=XLIM 

XG2=XGIT(NGIT) 

AROIP =ROGIT (NPO)+ROLIM5(XLIMXX3/3.DO 
AROUIP=ROUGIT(NPO) 

AEIP =EGIT (NPO)+ ELIM3(XLIMXX3/3.D0 

XP=X(2)/XCHARG 

DO 100 1=2, LL 

IP=I+1 

XI = XP 

AROI =AROIP 
AROUI=AROUIP 
AEI =AEIP 
XP=X(IP)/XCHARG 

IF(DABS(XP-1 .DO) .LT.l.D-10) XP=1 .DO 
CSOF=(XP.GE.l .DO) 

IF(XP.GE.XLIM) GO TO 101 
UNIFORM FLOW REGION 

DELVOL=(XLIM-XP)5<(XLIM5(3<2 + XLIM3«XP+XP3<3(2)/3.DO 
AROIP =ROGIT (NPO)+ROLIM5«DELVOL 
AROUIP=ROUGIT(NPO) 

AEIP =EGIT (NPO)+ ELIM3«DELV0L 
GO TO 102 
,01 CONTINUE 
NON UNIFORM FLOW REGION. 

IF( .NOT.CSOF) GO TO lOA 

LAST POINT. (THIS IS THE DETONATION FRONT POINT X=l). 
AROIP= 0. 

AROUIP=0. 

AEIP= 0. 

GO TO 102 
■OA CONTINUE 

IF(XP.LE.XG2) GO TO 103 
NGIT=NGIT-1 

IF(NGIT.LE.O) CALL SOFCBEGIN lOA. NGIT.LE.O.') 

XG1=XG2 

XG2=XGIT(NGIT) 

GO TO lOA 
.03 CONTINUE 

FRAC=(XP-XG1)/(XG2-XG1) 

IF(FRAC.LT.O. ) CALL SOFCBEGIN 103. FRAC.LT.O.') 
IF(FRAC.GT.l.DO) CALL SOFCBEGIN 103. FRAC.GT.l.') 

AROIP =(1 .DO-FRAC)5(ROGIT (NGIT+1 ) + FRAC5(R0GIT (NGIT) 
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CHA0299 
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CHA0312 

CHA0513 
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CHA0316 

CHA0317 

CHA0318 

CHA0319 

CHA0320 

CHA0321 

CHA0322 

CHA0323 

CHA032A 

CHA0325 

CHA0326 

CHA0327 

CHA0328 

CHA0329 

CHA0330 

CHA0331 

CHA0332 

CHA0333 

CHA033A 

CHA0335 

CHA0336 

CHA0337 

CHA0338 

CHA0339 

CHA03A0 

CHA03A1 

CHA03A2 

CHA03A3 

CHA03AA 

CHA0345 

CHA0346 

CHA0347 

CHA0348 

CHA0349 

CHA0350 

CHA0351 

CHA0352 

CHA0353 

CHA0354 

CHA0355 

CHA0356 

CHA0357 

CHA0353 

CHA0359 

CHA0360 

CHA0361 

CHA0362 

CHA0363 
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AR0UIP = (1 .DO-FRAC)3«ROUGIT(NGIT+1) + FRAC3(ROUGIT(NGIT) CHA0368 

AEIP =(1.D0-FRAC)5(EGIT (NGIT+1)+FRAC5«EGIT (NGIT) CHA0369 

102 CONTINUE CHA0370 

C COMPUTE MASS, MOMENTUM AND ENERGY DENSITIES. CHA0371 

IF(XP.LE.XLIM) GO TO 105 CHA0372 

C CONSERVATION-FORM DEFINITION OF MASS, MOMENTUM AND ENERGY DENSITY. CHA0373 
DVOL = (XP-XI)3«(XP3«5«2+XP5«XI+XI5«3«2)/3.DO CHA0374 

RO (I)=RONORM5t(AROI - AROIP)/DVOL CHA0375 

TENA(I)=RUNORM5((AROUI-AROUIP)/DVOL CHA0376 

E (I) = ENORM 5«(AEI - AEIPl/DVOL CHA0377 

GO TO 106 CHA0378 

105 CONTINUE CHA0379 

C UNIFORM FLOW REGION CHA0380 

RO (I) = RONORM5(ROLIM CHA0381 

TENA(I)=0. CHA0382 

E (I) = ENORM ELIM CHA0383 

106 CONTINUE CHA0384 

U(I)=TENA(I)/RO(I) CHA0385 

P(I) = (GAMA-1 .D0)5((Ea)-0.5D05(R0a)5(U(I)*x2) CHA0386 

PRINT lll,I,CSOF,U(I),Pa),ROa),E(I) CHA0387 

111 F0RMAT(/1X, *I,CS0F,U,P,R0,E=M4,L3,4D14.4) CHA0388 

IF(CSOF) GO TO 109 CHA0389 

100 CONTINUE CHA0390 

109 CONTINUE CHA0391 

DO 4 1=2, LL CHA0392 

Dxsia)=(xa+i)-x(i))5«RO(i) CHA0393 

4 CONTINUE CHA0394 

RETURN CHA0395 

end CHA0396 



SUBROUTINE TAYLOR(GAMA) 
IMPLICIT REAL5«8(A-H,0-Z,$) 



TAJ L O'R 



TAYLOR " SELF SIMILAR SPHERICAL DETONATION (CJ) FLOW FIELD 

COMMON /GGGG/G, G1 , G2 , G3 , G4 , G5 , G6 , G7 , G8 , G9 , G1 0 
COMMON /PAR/RHOO , QO , ROCJ , DCJ , UCJ , PCJ , DPSI , PSIMAX, CO . UO 
COMMON /GITN/NPO 

COMMON/GIT/ROLIM, EL IM,XGIT( 200 ), ROGITC 200 ), ROUGH (200 ), EGITC 200 ) 



CHA0397 

CHA0398 

CHA0399 

CHAOAOO 

CHAOAOl 

CHA0A02 

CHA0A03 

CHA040A 

CHA0405 






G=GAMA 
PRINT 101 
101 FORMAT(»l*) 

PRINT 110 

110 FORMATCIX, »G. I. TAYLOR SOLUTION. 
CALL INIDAT 
X=1.D0 
Y=0. 

U = U0 
C=C0 
AM=0. 

AT = 0 . 

AE=0. 

PSI=-DLOG(U) 

DO 1 N=1,NP0 



CHA0A07 

CHA0A08 

CHA0A09 

CHAOAIO 

N,PSI,U,C,X/AM,AT,AE= V/) CHAO All 

CHA0A12 
CHA0A13 
CHAOAIA 
CHA0A15 
CHA0A16 
CHA0A17 
CHA0A18 
CHA0419 
CHA0A20 
CHA0A21 
CHA0A22 
CHA0A23 
CHA0A2A 
CHA0A25 
CHA0426 
CHA0427 
CHA0428 
CHA0429 
CHA0430 
CHA0431 
CHA0432 
CHA0433 
CHA0434 
CHA0435 
CHA0436 
CHA0437 
CHA0438 
CHA0439 



XGIT (N)=X 
ROGIT (N)=AM 
ROUGIT(N)=AT 
EGIT (N)=AE 

PRINT 11, N,PSI,U,C,X,AM,AT,AE 

11 FORMAT (IX, I4,4D14.5/5X,3D14.5) 

CALL RUNGE(N,PSI,X,C,AM,AT,AE,PSIN,XN,CN,AMN,ATN,AEN) 

PSI=PSIN 

U=DEXP(-PSI) 

X = XN 

C=CN 

AM=AMN 

AT=ATN 

AE=AEN 

1 CONTINUE 

ROLIM=(C/CO)5^^G3 

ELIM=G55^(C/C0)5(XG4 

AMO =AM+( C/CO )5^^G35^X5^ 5^3/3. DO 
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AM0=AM0X3.D03«(G+1.D0)/G 

AE0=AE+(G5X(C/C0)3O«GA+0 .5D05<(C/C0)XXG33(U5<x2)xX5<x3/3.D0 
AE0=AE0X6 .D0X(G+1.D0)X(GXX2-1 .DO)/G 
PRINT 22,AM0,AE0 

FORMATC///1X, 'MASS AND ENERGY CHECK (SHOULD BE 1.)'// 

1 IX, »M0=',D17 .8,5X, *E0=’,D17 .3//) 

RETURN 
END 



CHA04A0 

CHAOAAl 

CHA0AA2 

CHA0AA3 

CHAOAAA 

CHA0AA5 

CHA0AA6 

CHA04A7 



SUBROUTINE INIDaT lUmr^ CHA0448 

IMPLICIT REAL5«8(A-H,0-Z,$) ' CHA0449 
COMMON /GGGG/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10 CHA0450 
COMMON /PAR/RHOO,QO,ROCJ,DCJ,UCJ,PCJ,DPSI,PSIMAX,CO,UO CHA0451 
COMMON /GITN/NPO CHA0452 



NP0=200 
PSIMAX=10 .DO 
UO=l .DO/(G+l .DO) 

CO=l.DO-UO 

DPSI=PSIMAX/DFLOAT(NPO) 

Gl=G-l.DO 

G2=G1/2.D0 

G3=2.D0/(G-1.D0) 

G4=2.D05(G/(G-1.D0) 

G5=G/((G+1 .D0)XX2X(G-1 .DO)) 

RETURN 

END 



J?ur(G£ 



SUBROUTINE RUNGE( N, PSI , X, C, AM, AT, AE, PSIN, XN, CN, AMN, ATN, AEN) 
IMPLICIT REAL5(8(A-H,0-Z,$) 

COMMON /GGGG/G, G1 , G2, G3 , G4, G5, G6 , G7 , G8 , G9 , G1 0 
COMMON /PAR/RHOO,QO,ROCJ,DCJ,UCJ,PCJ,DPSI,PSIMAX,CO,UO 
COMMON /GITN/NPO 



CHA0454 

CHA0455 

CHA0456 

CHA0457 

CHA0458 

CHA0459 

CHA0460 

CHA0461 

CHA0462 

CHA0463 

CHA0464 

CHAn44»i. 



CHA0466 

CHA0467 

CHA0468 

CHA0469 

CHA0470 






H=DPSI 
H2=H/2.D0 
H6=H/6 .DO 

CALL DERIV(PSI,X,C,AM,AT,AE, 

1 DXDP1,DCDP1,DMDP1,DTDP1,DEDP1) 

CALL DERIV(PSI+H2,X+H2^DXDP1,C+H23(DCDP1,AM,AT,AE, 
1 DXDP2,DCDP2,DMDP2,DTDP2,DEDP2) 

CALL DERIV(PSI+H2,X+H2^DXDP2,C+H23(DCDP2,AM,AT,AE, 
1 DXDP3,DCDP3,DMDP3,DTDP3,DEDP3) 

CALL DERIV(PSI + H,X+HXDXDP3,C+H5(DCDP3,AM,AT, AE, 

1 DXDP4,DCDP4,DMDP4,DTDP4,DEDP4) 

PSIN=PSI+H 

XN=X+H6X(DXDP1+2.D05((DXDP2+DXDP3) + DXDP4) 
CN=C+H63«(DCDPl + 2. D03<(DCDP2+DCDP3) + DCDP4) 

AMN = AM+H65«(DMDP1 + 2.D03«(DMDP2+DMDP3) + DMDP4) 

ATN = AT+H6K(DTDPl + 2. D03«(DTDP2+DTDP3) + DTDP4) 
AEN=AE+H65<( DEDPl + 2 . DOM( DEDP2+DEDP3) + DEDP4 ) 

RETURN 
END 






SUBROUTINE DERIV(PSI,X,C,AM,AT, AE, DXDP, DCDP, DMDP, DTDP, DEDP) 
IMPLICIT REALX8(A-H,0-Z,$) 

COMMON /GGGG/G, G1 , G2 , G3 , G4 , G5 , G6 , G7 , G8 , G9 , G1 0 
COMMON /PAR/RHOO,QO,ROCJ,DCJ,UCJ,PCJ, DP5I,PSIMAX,C0,U0 
COMMON /GITN/NPO 



CHA0472 

CHA0473 

CHA0474 

CHA0475 

CHA0476 

CHA0477 

CHA0478 

CHA0479 

CHA0480 

CHA0481 

CHA0432 

CHA0483 

CHA0484 

CHA0485 

CHA0486 

CHA0487 

CHA0488 

CHA0489 

CHA0490 



CHA0491 

CHA0492 

CHA0493 

CHA0494 

CHA0495 



U=DEXP(-PSI) CHA0497 

DXDP = 0 . 5D0^X5<(C-U+X)X(C+U-X)/D(5<2 CHA0498 

DCDP = -G25(UX(X-U)/C CHA0499 

DMDP = -(C/C0)XXG3XX?(X2?«DXDP CHA050 0 

DTDP=DMDP^U CHA0501 

DEDP = -( G55«( C/CO ) X3(G4+0 .SDOXC C/CO) ^XG3XU^X2)5<X^?(2^DXDP CHA0502 

RETURN CHA0503 

END CHAQ - S J)4- 



DOUBLE PRECISION FUNCTION RATIO(X) KATIO CHA0505 
IMPLICIT REAL)«8(A-H,0-Z,$) CHA0506 

RATIO=0. CHA0508 
IF(X.LE.1.D-8)RETURN CHA0509 
RATIO=2.DO/X CHA0510 
RETURN CHA0511 
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END 



CHA0512. 



DOUBLE PRECISION FUNCTION CROSS(X) CROSS CHA0513 

IMPLICIT REALX8(A-H,0-Z,$) CHA051A 

C CR0SS=1.D0 CHA0516 

CR0SS=XX3(2 CHA0517 

RETURN CHA0518 

END CHAO»>iq 



SUBROUTINE CYCEUL CMCet^L. CHA0520 

(L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN, CHA0521 

US,PS,UIDOT,PIDOT, CHA0522 

FIMZ,ZMDOT, CHA0523 

TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ) CHA0524 

IMPLICIT REALX8(A-H,0-Z,$) CHA0525 

DIMENSION X(L),U(L),P(L),RO(L),G(L),E(L),DU(L),DP(L),DRO(L), CHA0526 

DG(L),DXSI(L),MIN(L), CHA0527 

US(L),PS(L),UIDOT(L),PIDOT<L) CHA0528 

,TENA(L),FIRO(L),FIM(L),FIE(U CHA0529 

,GIP(L),VOL(U,Z(L),DZ(L) CHA0530 

,FIMZ(L),ZMDOT(L) CHA0S31 

COMMON /AB/AC50) CHA0532 

EQUIVALENCE ( LL, A( 2) ) , (T, A( 3) ) , (DT, A< A) ) , (COLELA, A( 13) ) CHA0S33 

EQUIVALENCE (KEYEK,A(16)) CHA0534 

EQUIVALENCE ( STAB, A( 18 ) ) , ( DTBA, A( 19) ) , ( DTKOD, A( 20 ) ) , (KDT, A(21) ) CHA0S35 
COMMON /GAM/GAMA, NG,MU2,G1,G2,G3,GA,G5,G6,G7,G8,G9,G10,G11 CHA0536 

,G12,G13,G1A,G15,G16,G17,G18,G19,G20,G21,G22,G23 CHA0537 

,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G3A,G35 CHA0538 

REALX8 NG,MU2 CHA0539 

COMMON /TOT/AMTOT,ETOT,EKTOT,EPTOT,TENTOT CHA05A0 

COMMON /AZOV/ISAFA,NORIMN,USAF,PSAF,ROSAF,GSAF,ESAF,DPSAF CHA05A1 

,DXSIL,DXSIR CHA05A2 

LOGICAL NORIMN CHA05A3 

COMMON /STEPO/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR, CHA054A 

RSTARL,RSTARR,GSTARL,GSTARR, CHA0595 

CL,CR,CSTARL,CSTARR,SL,SR,NL,WR,UW(6) CHA0546 

,LAMDAL,LAMDAR,RATEL, RATER, TEMPL,TEMPR,TEMPSL,TEMPSR CHA0547 
,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR CHA0543 

REALX8 LAMDAL,LAMDAR CHA0549 

LOGICAL HELEML,HELEMR CHA0550 

COMMON /STEP1/DUIDT,DPIDT,DGIDTL, DGIDTR,DRIDTL,DRIDTR CHA0551 

,ASTARL,ASTARR,LAMDSL, LAMDSR,DSDAL,DSDAR,DZDAL,DZDAR CHA0552 

,RAT,SH CHA0553 

, BETACL , BETACR, DSDASL , DSDASR , DZDASL , DZDASR CHA0554 

REALX8 LAMDSL,LAMDSR,DSDAL,DSDAR,DZDAL,DZDAR CHA0555 

COMMON /GRADS/DUDXIL, DPDXIL, DGDXIL, DRDXIL, DZDXIL, DSDXIL, CHA0556 

1 DUDXIR,DPDXIR,DGDXIR,DRDXIR,DZDXIR,DSDXIR CHA0557 

COMMON /FI/FIH1,FIH2,FIH3,UXN,PXN,GXN,R0XN,ZXN CHA0558 

1 ,GIH CHA0559 

2 ,FIH4,ZMD0TL,ZMD0TR CHA0560 

COMMON/DETO/QDET,PCJDET,RCJDET,UCJDET,DCJDET,PODET,ROODET, CHA0561 

1 RATE,TEMPC CHA0562 

C0MM0N/PULS/PRESSa0),PULSEl(10),PULSE2(10),PULSE3(10) ,PULSE4(10) CHA056 3 
DATA ERRP/O.DO/ CHA0564 

DATA NERRP/0/ CHA0565 

C DATA KOTZ/7777777777B/ CHA0566 



DT2=DT/2.D0 


CHA0568 


UXN=0. 


CHA0573 


PXN=0. 


CHA0574 


R0XN=0. 


CHA0575 


ZXN=0. 


CHA0576 


DO 1 1=2, L 


CHA0577 


IM=I-1 


CHA0578 


UXNM=UXN 


CHA0579 


PXNM=PXN 


CHA0580 


ROXNM=ROXN 


CHA0581 


ZXNM=ZXN 


CHA0582 


UL=U(IM)+0 .SDOXDUCIM) 


CHA0583 


PL=P(IM)+0.5D0XDP(IM) 


CHA0584 


ROL=RO(IM) + 0 . SDOXDROdM) 


CHA0585 


GL=DSQRT(GAMAXPLXR0L) 


CHA0586 


CL=GL/ROL 


CHA0587 
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2L=Z(IM)+0.5D0XDZ(IM) 


CHA0588 


IF(ZL .GT.l .DO) ZL=l.DO 


CHA0589 


IF(ZL.LT.O. ) ZL=0. 


CHA0590 


SL=PL/(G1XR0LXXGAMA) 


CHA0591 


TEMPL=PL/ROL 


CHA0592 


UR=U(I)-0.5D0XDU(I) 


CHA0593 


PR=P(I)-0.5D0XDP(I) 


CHA059A 


ROR=RO(I)-0.5DOXDRO(I) 


CHA0595 


GR=DSQRT(GAMAXPRXR0R) 


CHA0596 


cr=gr/ror 


CHA0597 


ZR=Z(I)-0.5D0XDZ(I) 


CHA0598 


IF(ZR.GT.l.DO) ZR=l.DO 


CHA0599 


IF(ZR.LT.O. ) ZR=0. 


CHA0600 


SR=PR/(G1KR0RXXGAMA) 


CHA0601 


TEMPR=PR/ROR 


CHA0602 

CHA0603 


CALL RIEMAN(L,I,MIN) 


CHA0604 

CHA0605 


DUDXIL=DU(IM)/DXSI(IM) 


CHA0606 


DPDXIL=DP(IM)/DXSI(IM) 


CHA0607 


DRDXIL=DRO(IM)/DXSI(IM) 


CHA0608 


DGDXIL=0.5D0XGLX(DPDXIL/PL+DRDXIL/R0L) 


CHA0609 


DZDXIL=DZ(IM)/DXSI(IM) 


CHA0610 


DSDXIL=SLX(DPDXIL/PL-GAMAXDRDXIL/ROL) 


CHA0611 


DUDXIR=DU(I)/DXSI(I) 


CHA0612 


DPDXIR=DP(I)/DXSI(I) 


CHA0613 


DRDXIR=DRO(I)/DXSI(I) 


CHA0614 


DGDXIR=0 . 5D0XGRX( DPDXIR/PR+DRDXIR/ROR) 


CHA0615 


DZDXIR=DZ(I)/DXSI(I) 


CHA0616 


DSDXIR=SRx(DPDXIR/PR-GAMAXDRDXIR/ROR) 


CHA0617 


SH=CROSS(X(D) 


CHA0618 


RAT=RATIO(X(D) 


CHA0619 

CHA0620 


CALL MAGA(L,I,MIN) 


CHA0621 

CHA0622 


US(I)=USTAR 


CHA0623 


PS(I)=PSTAR 


CHA0624 


UIDOT(I)=DUIDT 


CHA0625 


PIDOT(I)=DPIDT 


CHA0626 

CHA0627 


CALL FLUXE(L,I,MIN) 


CHA0623 

CHA0629 


FIR0(I)=FIH1 


CHA0630 


FIM (I)=FIH2 


CHA0631 


FIE (I)=FIH3 


CHA0632 


FIMZ(I)=FIH4 


CHA0633 


GIP(I)=GIH 


CHA0634 


DU(IM)=UXN-UXNM 


CHA0635 


DP(IM)=PXN-PXNM 


CHA0636 


DRO(IM)=ROXN-ROXNM 


CHA0637 


DZ(IM)=ZXN-ZXNM 


CHA0638 


STATIONS OUTPUT 


CHA0639 


IF((I-42)3«(I-62)x(I-82)3f(I-102) .NE.O) GO TO 1 


CHA0640 


NPU = 0 


CHA06A1 


IF(I.EQ.42) NPU=1 


CHA0642 


IF(I.EQ.62) NPU=2 


CHA06A3 


IF(I.EQ.82) NPU=3 


CHA06A4 


IF(I.EQ.102)NPU=4 


CHA0645 


IF(NPU.EQ.O) CALL SOFCFLUXE 90. NPU.EQ.O') 


CHA0646 


PRESS(NPU)=GIH+FIH2 


CHA0647 


PULSEI(NPU)=PULSE1(NPU) + DT5«GIH 


CHA0643 


PULSE2(NPU)=PULSE2(NPU) + DT5t(GIH+FIH2) 


CHA0649 


PULSE3(NPU)=PULSE3(.NPU) + DT5(FIH1^CR0SS(X(I) ) 


CHA0650 


PULSE4(NPU)=PULSE4(NPU) + DT3fFIH25«CROS3(X(I) ) 


CHA0651 


CONTINUE 


CHA0652 

CHA0653 


AMTOT=0 . 


CHA0654 


ETOT=0. 


CHA0655 


EKTOT=0. 


CHA0656 


EPTOT=0. 


CHA0657 


TENTOT=0. 


CHA0653 


FI1=FIR0(2) 


CHA0659 



37 



FILE; CHARGER FORTRAN 



A1 



C 

C 

291 

C 

29 

2 

C 

200 



FI2=FIM (2) 

FI3=FIE (2) 

FIA=FIMZ(2) 

GI2=GIP(2) 

SH=CR0SS(X(2)) 

DO 2 1=2, LL 

IP=I+1 

FIM1=FI1 

FIM2=FI2 

FIM3=FI3 

FIMA=FI4 

GIM2=GI2 

SHM=SH 

FI1=FIR0(IP) 

FI2=FIM (IP) 

FI3=FIE (IP) 

FI4=FIMZ(IP) 

GI2=GIP (IP) 

SH=CROSS(X(IP)) 

DVOL=VOL(I) 

ROOLD=RO(I) 

POLD=P(I) 

EOLD=E(I) 

UOLD=U(I) 

ZOLD=Z(I) 

ZK0DM=Z0LDXR00LD 

TOLD=POLD/ROOLD 

DX=X(IP)-X(I) 

DTVOL=DT/DVOL 

RO(I)=RO(I)-DTVOL5((SH5(FIl-SHM3(FIMl) 

TENA(I)=TENA(I)-DTV0LX(SH3(FI2-SHMXFIM2)-(DT/DX)X(GI2-GIM2) 
E(I) = E(I)-DTVOL3«(SH5(FI3-SHM5(FIM3) 

U(I)=TENA(I)/RO(I) 

Z(I) = (ZKODM-DTVOL5«(SHXFIA-SHM5(FIMA))/RO(I) 

IF(Z(I) .GT.l.DO) Z(I)=1.D0 
IF(Z(I).LT.O.) Z(I)=0. 

UAV=U(I) 

ROAV=RO(I) 

EP = E( I )-0 . 5D05«ROAV3(UAV3(3«2 

IF(EP.GT.O.) GO TO 291 

NERRP=NERRP+1 

ERRP = ERRP+(1 .D-3-EP)5«DVOL 

IF(ERRP.GT.0.2AD0) GO TO 291 

EP=l.D-a 

CONTINUE 

IF(EP.LE.O.) GO TO 7001 
P(I)=G15«EP 

G(I) = DSQRT(GAMA3fP(I)XR0(I)) 

UPC=DABS(U(I))+G(I)/RO(I) 

DTI=STAB3«DX/UPC 

IF(DTI .GT.DTBA) GO TO 29 

DTBA=DTI 

KDT = I 

CONTINUE 

DXSI(I)=R0(I)5«DX 

ET0T = ET0T+E(I)3«DV0L 

EPT0T = EPT0T+EP3<DV0L 

AMTOT = AMTOT+RO(I)>«DVOL 

TENTOT=TENTOT+TENA(I)XDVOL 

CONTINUE 

EKTOT=ETOT-EPTOT 

IF(COLELA. EQ. 0. ) GO TO 200 
CALL DCOLE(L,X,U ,DU ,MIN,1) 

CALL DC0LE(L,X,P,DP,MIN,2) 

CALL DC0LE(L,X,R0,DR0,MIN,3) 

CALL DCOLE(L,X,Z,DZ,MIN, A) 

CONTINUE 

CALL BD0K1(L,X,U , DU ,MIN,1) 



CHA0660 

CHA0661 

CHA0662 

CHA0663 

CHA066A 

CHA0665 

CHA0666 

CHA0667 

CHA0668 

CHA0669 

CHA0670 

CHA0671 

CHA0672 

CHA0673 

CHA067A 

CHA0675 

CHA0676 

CHA0677 

CHA0678 

CHA0679 

CHA0680 

CHA0681 

CHA0682 

CHA0683 

CHA068A 

CHA0685 

CHA0686 

CHA0687 

CHA0688 

CHA0689 

CHA0690 

CHA0691 

CHA0692 

CHA0693 

CHA069A 

CHA0695 

CHA0696 

CHA0697 

CHA0698 

CHA0699 

CHA0700 

CHA0701 

CHA0702 

CHA0703 

CHA070<4 

CHA0705 

CHA0706 

CHA0707 

CHA0708 

CHA0709 

CHA0710 

CHA0711 

CHA0712 

CHA0713 

CHA071A 

CHA0715 

CHA0716 

CHA0717 

CHA0718 

CHA0719 

CHA0720 

CHA0721 

CHA0722 

CHA0723 

CHA072A 

CHA0725 

CHA0726 

CHA0727 

CHA0728 

CHA0729 

CHA0730 

CHA0731 
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CALL BD0K1(L,X,P,DP,MIN,2) 

CALL BD0K1(L,X,R0,DR0,MIN,3) 

CALL BDOKl(L,X,Z,DZ,MIN,A) 

PRINT 901,(NN,PRESS(NN),PULSE1(NN),PULSE2(NN), 

1 PULSE3(NN)/AMT0T,PULSEA(NN)/TENT0T,NN=1,4) 

101 F0RMAT(1X,2( , 13, 5D1 1 . 3, '/')/) 

IF(DABS(T-A(5)).LT.l.D-6) PRINT 91 1 , NERRP, ERRP 
111 FORMATC//1X, 'NERRP,ERRP=M5,D15.5/) 

RETURN 

001 CONTINUE 

PRINT 7101, I,ROAV,UAV,DRO(I),DU(I),E(I),EP,ZNEW,ZNEW-1.DO,EPI 
101 FORMATC//1X, 'FROM CYCEUL . NEGATIVE EP . IN CELL I=',I6// 

1 IX, 'ROAV,UAV,DRO(I),DU(I)=',AD18.8// 

2 IX, 'E(I),EP,ZNEW,ZNEW-1,EPI=',5D14.6//) 

CALL PRINT 



(L,X,U,P,RO,G, E, DU,DP,DRO,DG,DXSI,MIN, 
US,PS,UIDOT,PIDOT, 

FIMZ,ZMDOT, 

TENA,FIRO, FIM,FIE,GIP,VOL,Z,DZ) 

CALL SOF( 'CYCEUL 7001, NEGATIVE EP') 

RETURN 

END 

SUBROUTINE' SAFAE 



1 (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN, 

2 US,PS,UIDOT,PIDOT, 

X FIMZ,ZMDOT, 

3 TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ) 

IMPLICIT REAL5«8(A-H,0-Z,$) 

DIMENSION X(L),U(L),P(L),RO(L),G(L),E(L),DU(L),DP(L),DRO(L), 

1 DG(L),DXSI(L),MIN(L), 

2 US(L),PS(L),UIDOT(L),PIDOT(L) 

3 ,TENA(L),FIRO(L),FIM(L),FIE(L) 

4 ,GIP(L),VOL(L),Z(L),DZ(L) 

5 ,FIMZ(L),ZMDOT(L) 

COMMON /AB/A(50) 

EQUIVALENCE ( LL , A(2) ) , ( T, A( 3 ) ) , ( DT, A( 4) ) , ( NCYC, A( 12 ) ) 
EQUIVALENCE (UGAL,A(15)) 

COMMON /GAM/GAMA, NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11 

1 ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23 

2 ,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,G35 
REALMS NG,MU2 

COMMON/DETO/QDET,PCJDET,RCJDET,UCJDET,DCJDET,PODET,ROODET, 

1 RATE,TEMPC 

C0MM0N/DIFFUS/U2,P2,R02,ARW 



CHA0732 

CHA0733 

CHA0734 

CHA0735 

CHA0736 

CHA0737 

CHA0738 

CHA0739 

CHA0740 

CHA0741 

CHA0742 

CHA0743 

CHA0744 

CHA0745 

CHA0746 

CHA0747 

CHA0748 

CHA0749 

CHA0750 

CHA0751 

CHA0752 

CHA0753 

CHA0754 

CHA0755 

CHA0756 

CHA0757 

CHA0758 

CHA0759 

CHA0760 

CHA0761 

CHA0762 

CHA0763 

CHA0764 

CHA0765 

CHA0766 

CHA0767 

CHA0768 

CHA0769 

CHA0770 

CHA0771 

CHA0772 

CHA0773 

CHA0774 

CHA0775 






RIGID B.C. AT 1=2 

UCn=-U(2) 

PCn=P(2) 

G(1)=G(2) 

R0(1)=R0(2) 

Z(n=Z( 2 ) 

DU(1)=DU(2) 

DP(1)=-DP(2) 

DG(1)=-DG(2) 

DR0(1)=-DR0(2) 

DXSI(1)=DXSI(2) 

OUTFLOW B.C. AT I=L 

U(U = U(LU + DU(LL)/2 . DO 
P(U=P(LU + DP(LU/2,D0 
ROCL)=ROCLU + DRO(LL)/2.DO 
G(U=G(LU + DG(LU/2. DO 
Z(U=Z(LU + DZ(LU/2. DO 
DU(U=0. 

DPCU = 0. 

DGCL)=0. 

DRO(U = 0. 

DZ(U=0 . 

DXSI(U=DXSI(LL) 
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RETURN 

END 
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SUBROUTINE BDOKl(L,X, V,DV,MIN,NV) 

IMPLICIT REALX8(A-H,0-2,$) 

DIMENSION X(L),V(U,DV(L),MIN(L) 

COMMON /AB/A(50) 

EQUIVALENCE ( LL , A( 2) ) , ( KEYMON, A( 1 1 )) 

COMMON /DRAW/GODELX, GODELY, UMIN, UMAX, PMIN, PMAX, ROMIN, ROMAX 
1 ,XMIN,XMAX,SMIN,SMAX,IVERSA 

COMMON /M0NIT/CASEAV(4),NC1A(4),NF16(6), 

1 NMONU(A),NMONP(A),NMONRO(A),NMONZ(A) 

DIMENSION NM0NV(A,4) 

EQUIVALENCE (NMONVd, l),NMONU(l) ) 

DIMENSION NAMEV(A) 

DATA NAMEV/'U', 'P', 'RO', *ZV 
DATA EPS/1, D-9/ 
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GO TO (1,2, 3, 4), NV 


CHA0822 


1 


AMIDA=(UMAX-UMIN)XX2 


CHA0823 




GO TO 9 


CHA0824 


2 


AMIDA=(PMAX-PMIN)XX2 


CHA0825 




GO TO 9 


CHA0826 


3 


AMIDA=(R0MAX-R0MIN)XX2 


CHA0827 




GO TO 9 


CHA0828 


4 


AMIDA=1.D0 


CHA0829 




GO TO 9 


CHA0830 


9 


CONTINUE 


CHA0831 




AMIDA=AMIDAXEPSXX2 


CHA0832 




EPSA=DSQRT(AMIDA) 


CHA0833 




DO 29 1=2, LL 


CHA0834 




ICAT=0 


CHA0835 




IF(DABS(DV(D) .LE.EPSA) DV(I) = 0. 


CHA0836 




IF(DV(I).EQ.O.) GO TO 29 


CHA0837 




VLEFT=V(I)-0.5D05(DV(I) 


CHA0838 




VRIGHT=V(I)+0.5D0XDV(I) 


CHA0839 




VM=V(I-1) 


CHA0840 




VP=V(I+1) 


CHA0841 




SIGN=(VP-V(I))X(V(I)-VM) 


CHA0842 




IF(SIGN.GT.-AMIDA) GO TO 22 


CHA0843 


21 


DV(I)=0. 


CHA0844 




ICAT=1 


CHA0845 




GO TO 20 


CHA0846 


22 


CONTINUE 


CHA0847 




SIGN=(VP-VM)XDV(I) 


CHA0848 




IF(SIGN.GT.-AMIDA) GO TO 24 


CHA0849 


23 


DV(I)=0.5D0*(VP-VM) 


CHA0850 




VLEFT=V(I)-0.5D0XDV(I) 


CHA0851 




VRIGHT=V(I)+0.5D0XDV(I) 


CHA0852 




ICAT=2 


CHA0853 


24 


SIGN=(VLEFT-VM)3(DV(I) 


CHA0854 




IF(SIGN.GT.-AMIDA) GO TO 26 


CHA0855 


25 


VLEFT=VM 


CHA0856 




VRIGHT=2.D0*V(I)-VLEFT 


CHA0857 




DV(I)=VRIGHT-VLEFT 


CHA0858 




ICAT=3 


CHA0859 


26 


SIGN=(VP-VRIGHT)5(DV(I) 


CHA0860 




IF (SIGN.GT.-AMIDA) GO TO 28 


CHA0861 


27 


VRIGHT=VP 


CHA0862 




VLEFT=2.D0XV(I)-VRIGHT 


CHA0863 




DV(I)=VRIGHT-VLEFT 


CHA0864 




ICAT=3 


CHA0865 


28 


IF(DABS(DV(D) .LE.0.5D05CDABS(VP-VM)) GO TO 31 
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30 


DV(I)=0.5D0X(VP-VM) 
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ICAT=4 


CHA0868 


31 


CONTINUE 


CHA0869 


20 


CONTINUE 


CHA0870 




IF (DABS(DV(I)).GT.EPSA) GO TO 40 


CHA0871 




DV(I)=0. 


CHA0872 


40 


CONTINUE 


CHA0873 




IF (ICAT.GT.O) NMONV(ICAT,NV)=NMONV(ICAT,NV)+1 
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29 


CONTINUE 
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RETURN 

END 
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SUBROUTINE DCOLE( L , X, V, DV, MIN, NV) 
IMPLICIT REALX8(A-H,0-Z,$) 
DIMENSION X(L),V(L),DV(L),MIN(L) 
COMMON /AB/A(50) 

EQUIVALENCE (LL,A(2)) 



J>COL£ 
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DO 1 1=2, LL 

IM=I-1 

IP=I+1 

DV(I)=0.5D0X(VaP)-VaM)) 
CONTINUE 
RETURN 
END 



SUBROUTINE PRINT 

(L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN, 
US,PS,UIDOT,PIDOT, 

FIMZ,ZMDOT, 

TENA,FIRO, FIM,FIE,GIP,VOL,Z,DZ) 

IMPLICIT REALX8(A-H,0-Z,$) 

DIMENSION X(L),U(L),P(L),RO(L),G(L),E(L),DU(L),DP(L),DRO(L), 
DG(L),DXSI(L),MIN(L), 

US(L),PS(L),UIDOT(L),PIDOT(L) 

,TENA(L),FIRO(L),FIM(L),FIE(L) 

,GIP(L),VOL(L),Z(L),DZ(L) 

,FIMZ(L),ZMDOT(L) 

/TOT/AMTOT,ETOT,EKTOT,EPTOT,TENTOT 
/STEPO/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR, 
RSTARL,RSTARR,GSTARL,GSTARR, 
CL,CR,CSTARL,CSTARR,SL,SR,WL,WR,UW(6) 
,LAMDAL,LAMDAR,RATEL, RATER, TEMPL,TEMPR,TEMPSL,TEMPSR 
,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR 
REALX8 LAMDAL,LAMDAR 
LOGICAL HELEML,HELEMR 
COMMON /AB/A(50) 

EQUIVALENCE (LL, A(2) ) , (T, A(3) ) , (NCYC, A(12) ), (DT,A(A) ) 

EQUIVALENCE (UGAL,A(15)) 

C0MM0N/DIFFUS/U2, P2, R02, ARM 

COMMON/DETO/QDET, PCJDET, RCJDET, UCJDET, DCJDET , PODET, ROODET , 

1 RATE,TEMPC 

COMMON /GAM/GAMA, NG, MUZ, G1,G2,G3,GA,G5,G6,G7,G8,G9,G10, Gil 

1 ,G12,G13,G1‘4,G15,G16,G17,G1S,G19,G20,GZ1,G22,G23 

2 , G2A , G25, G26 , G27 , G28 , G29 , G30 , G31 , G32, G33, G3A, G35 
REAL5(8 NG,MU2 

COMMON /M0NIT/CASEAV(A),NC1‘4(A),NF16(6), 

1 NMONUC A),NMONP(A),NMONRO(A),NMONZ(A) 

DIMENSION CASAVKA) 

LOGICAL FULLPR 
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FR N'T 



COMMON 

COMMON 



FULLPR=.TRUE. 

PRINT 1 
FORMAT(lHl) 

PRINT 2, T,DT,NCYC 

FORMAT ( IX, 1 OX, 'RESULTS AT T= ' , D1 1 . 5 , 5X, ' DT= ' , D1 1 . 5 , 5X, 'NCYC=> 
1 15//) 

PRINT 3, AMTOT,ETOT,EKTOT,EPTOT,TENTOT 

FORMAT ( IX, ' AMTOT= ' , D20 . 1 A , 2X, ' ETOT, EKTOT, EPTOT= ' , 3D22 . 1 A/ 

1 IX, 'TENTOT=' ,D21 .lA//) 



FORMATCIX, 



1 
2 
3 

AA FORMATCIX,' 
1 
2 
3 



X 

RO 

DU 

DG 

ZMDOT 

AMDOTN 

AMACH 



U 

G 

DP 

DZ' ) 
US 

FIMZ 

TEMP 

ENTRO 



I 

I 

j 

I I 

9 

I 

I 

I I 

, 

' ) 



DRO 



P 

Z 



PS 

AMDOT 

ENTALP 



FORMATCIX) 

IF CUGAL.NE.O,) PRINT 6, UGAL 
FORMATC/llX, 'INITIAL VELOCITY 
DO 10 1=1, L 

IF CM0DCI,10) ,NE.l) GO TO 11 



CORRESPONDS TO UGAL = ' , D15 . 6/ ) 
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PRINT 5 
PRINT 4 
PRINT 44 
PRINT 5 

11 CONTINUE 

PRINT 12,I,Xa),U(I),P(I),R0(I),G(I),Z(I),DU(I),DP(I),DR0a), 

1 DG(I),DZ(I) 

12 FORMATCIX, 13, 6012.5,5011.4) 

ENTR0=P(I)/R0(I)XXGAMA 

IF( .NOT.FULLPR) GO TO 131 
IF(I.EQ.l) GO TO 131 
IM=I-1 

UL=U(IM)+0 .5X0U(IM) 

PL=P(IM)+0.53«OP(IM) 

ROL=RO(IM)+0 .5XDR0CIM) 

GL=G(IM)+0.5xDG(IM) 

CL=GL/ROL 

ZL=Z(IM)+0.5XOZ(IM) 

IF(ZL.LT.O.) ZL=0. 

UR=U(I)-0.5XOU(I) 

PR=P(I)-0.5XOP(I) 

GR=G(I)-0.5XDG(I) 

ROR=ROa)-0.5XDRO(I) 

CR=GR/ROR 
ZR=Z(I)-0.5XDZa) 
ifczr.lt. 0.) ZR=0. 

IFCPL.LE.O.) PL=1.0-8 
IFCPR.LE.O.) PR=1.0-8 
CALL RIEMAN(L,I,MIN) 

XI=XCI) 

RSTAR=RSTARL 

IFCUSTAR.lt. 0.) RSTAR=RSTARR 
ZSTAR=ZL 

IFCUSTAR.LT.O. ) ZSTAR=ZR 
AMACH=USTAR/0SQRTCGAMAXPSTAR/RSTAR) 

AMD0T=RSTAR3<USTARXCR0SSCXI) 

IFCI.NE.2) GO TO 132 
AMOOTO=AMOOT 

I FC OABSC AMOOTO ) .LT. 1.0-12) AMOOTO=1 .00 
132 CONTINUE 

AMDOTN=AMDOT/AMOOTO 

ENTALP=CGAMA/CGAMA-1.00))XPSTAR/RSTAR+0 .500XUSTARXX2+QOETXZSTAR 
ARW=1 .00 

TEMP=PSTAR/CRSTARXARW) 

PRINT 13,USCI),PSCI), 

1 ZMOOTCI ) , FIMZCI) , AMOOT, AMOOTN, TEMP, ENTALP,AMACH, ENTRO 

13 F0RMATC4X,12X,5012.5,6011 .4) 

131 CONTINUE 

10 CONTINUE 
C JOB STATISTICS 
00 40 1=1,4 
CASAV1CI)=0. 

IF CNC14CI) .NE.O) CASAVl C I ) =CASEAVC I )/OFLOATC NC14C I ) ) 

40 CONTINUE 
PRINT 30 

30 F0RMATC///1X, IOC 'X* ),3X, ' JOB STATISTICS ', 3X, IOC ' X ')// ) 

PRINT 31, CNC14CI),I=1,4) 

31 FORMATCIX, ’NO. OF VARIOUS CASES IN RIEMAN SOLVER NC14CNCASE)= ' , 

1 4110) 

PRINT 301, CCASAV1CI),I=1,4) 

301 FORMATC/IX, 'AVERAGE NUMBER OF ITERATIONS IN RIEMAN SOLVER', 

1 IX,' CASAV1CNCASE)=',4CF6.2,4X)) 

PRINT 32, CNF16CI),I=1,6) 

32 FORMATC/IX, 'NO. OF VARIOUS FLUX CASES NF16CNFLUX)= ' ,6110) 
ICAT0=4 

PRINT 33, CNMONUCI),I = 1,ICATO),CNMONPCI),I = 1,ICATO), 

1 CNMONROCI), I=1,ICAT0), CNMONZCI), 1=1, ICATO) 

33 FORMATC/IX, 'NO. OF MONOTONICITY INTERVENTIONS FOR EACH VAR.', 

1 IX, 'IN EACH CATEGORY.'/ 

1 1X,'NM0NU CICAT)=' ,4110/ 

1 1X,'NM0NP CICAT)=' ,4110/ 

1 IX, 'NMONROCICAT)=',4I10/ 
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RETURN 
END 



1X,'NM0NZ (ICAT)=‘,AI10/) 



SUBROUTINE SOF(ISTOP) 

IMPLICIT REALX8(A-H,0-Z,$) 

DIMENSION ISTOP(l) 

PRINT 1, ISTOP 

F0RMAT(//1X, SHXXX, 3X, 20A4, 3X, 3HXXX///) 
PRINT 1 
XX=-1 .DO 
YY=DSQRT(XX) 

STOP 

DiP- 
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SUBROUTINE RIEMANC L , I , MIN) 

IMPLICIT REAL5(8(A-H,0-Z,$) 

DIMENSION MIN(L) 

COMMON /STEPO/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR, 

1 RSTARL,RSTARR,GSTARL,GSTARR, 

2 CL,CR,CSTARL,CSTARR,SL,SR,WL,WR,UW(6) 



«\sman 
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3 

A ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR CHA1317 

REALX8 LAMDAL,LAMDAR CHA1318 

LOGICAL HELEML,HELEMR CHA1319 

COMMON /STEP1/DUIDT,DPIDT,DGIDTL,DGIDTR,DRIDTL,DRIDTR CHA1320 

2 ,ASTARL,ASTARR,LAMDSL,LAMDSR,DSDAL,DSDAR,DZDAL,DZDAR CHA1321 

3 ,RAT,SH CHA1322 

A ,BETACL,BETACR,DSDASL,DSDASR,DZDASL,DZDASR CHA1323 

REALX8 LAMDSL,LAMDSR,DSDAL,DSDAR,DZDAL,DZDAR CHA132A 

COMMON /DRAW/GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX CHA1325 

1 ,XMIN,XMAX,SMIN,SMAX, IVERSA CHA1326 

COMMON /GAM/GAMA , NG, MUZ , G1 , G2 , G3 , GA , G5, G6 , G7 , G8 , G9 , G1 0 , G1 1 CHAl 327 

1 ,G12,G13,G1A,G15,G16,G17,G18,G19,G20,G21,G22,G23 CHA1328 

2 ,G2A,G25,G26,G27,G28,G29,G30,G31,G32,G33,G3A,G35 CHA1329 

REALX8 NG,MU2 CHA1330 

COMMON /AB/A(50) CHA1331 

COMMON /M0NIT/CASEAV(A),NC1A(A),NF16(6), CHA1332 

1 NMONU(A),NMONP(A),NMONRO(A),NMONZ(A) CHA1533 

DATA NMAX/63/ CHA1335 

DATA EPS/1. D-8/ CHA1336 

DATA NTRY/0/ CHA1337 






UW(6)=1 .D20 
WL=0. 

WR=0. 

ZETAL=PL3«XG8 

ZETAR = PR5«5«G8 

CLG=CL/GAMA 

CRG=CR/GAMA 

ZSTARL=ZL 

ZSTARR=ZR 

IF (ZETAL.lt. ZETAR) GO TO 102 
LEFT PRESSURE IS HIGHER 
)1 CONTINUE 

EVERR=(PL-PR)/PR 

USR=UR+CRGXEVERR/DSQRT( 1 . D0+G6XEVERR) 
SRR-USR 

UEL=UL-G75(CL5((ZETAR-ZETAU/ZETAL 

SLL=UEL 

NL=2 

NR = 2 

IF (USR.GE.UL) NL=1 
IF (UEL.LE.UR) NR=1 
IF (DABS(EVERR) .LT.EPS) GO TO 100 
IF (NL.EQ.2.AND.NR.EQ.1) GO TO 7001 
GO TO 100 

RIGHT PRESSURE IS HIGHER 
32 CONTINUE 

EVERL=(PR-PL)/PL 

USL = UL-CLG^EVERL/DSQRT(1 . D0+G65^EVERU 
SLL=USL 

UER = UR+G75(CRX(ZETAL-ZETAR)/ZETAR 
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SRR=UER 

NL=2 

NR=2 

IF (UER.GE.UL) NL=1 
IF (USL.LE.UR) NR=1 
IF (DABS(EVERL) .LT.EPS) GO TO 100 
IF (NL.EQ.1.AND.NR.EQ.2) GO TO 7001 
GO TO 100 
100 CONTINUE 

IF (NL .EQ.1.AND.NR.EQ.2) NCASE=1 
IF (NL.EQ.2.AND.NR.EQ.2) NCASE=2 
IF (NL.EQ.2.AND.NR.EQ.1) NCASE=3 
IF (NL.EQ.l.AND.NR.EQ.l) NCASE=A 

IF(DABS(PL-PR)+DABS(UL-UR) .LT.EPSX(PMAX-UMIN)) NCASE=4 
UMIDA=EPSXDMAX1 (CL , CR) 

DUDZL=-G7XCL/ZETAL 
DUDZR= G7XCR/ZETAR 

ZETA= ( -( UR-UL)+ZETAR3(DUDZR-ZETAL5(DUDZL )/ ( DUDZR-DUDZL ) 

IF (ZETA.LE.O.) GO TO 7002 

N=0 

GO TO (1,2, 3, 4), NCASE 
C THE CASE ES 

1 ITYPE=NCASE 
HELEML=. FALSE. 

HELEMR=.TRUE. 

11 N=N+1 

IF (N.GT.NMAX) GO TO 7003 
ZETAF=ZETA 

UEL=UL-G73(CL5«(ZETAF-ZETAL)/ZETAL 

PPR=(ZETAF/ZETAR)5(5(NG 

EVERR=PPR-1.D0 

SQRR=DSQRT(1.D0+G6XEVERR) 

USR=UR+CRG*EVERR/SQRR 

DU=UEL-USR 

IF (DABS(DU) .LE.UMIDA) GO TO 10 

DUDZR=NGXCRG5((PPR/ZETAF)5((1 .D0+G9XEVERR)/SQRRXX3 
ZETA=ZETAF+DU/ ( DUDZR-DUDZL ) 

GO TO 11 
10 CONTINUE 

USTAR=(UEL+USR)/2.D0 
IF(DABS(USTAR) .LT.EPSXUMAX) USTAR=0. 

PSTAR = PPRJ«PR 
CSTARL=CL+(UL-USTAR)/G7 
RSTARL=GAMAJ«PSTAR/CSTARLXX2 
GSTARL=CSTARL5fRSTARL 

C EQU. NO. 69.01 OF THE BOOK BY COURANT-FRIEDRICHS . 
WHR=G113((USTAR-UR)XR0R 
WR=WNR+DSQRT(GRXX2+WHRXX2) 

RSTARR = RORXWR/(WR-ROR5((USTAR-UR)) 

GSTARR = DSQRT(GAMA5«PSTAR5«RSTARR) 

CSTARR=GSTARR/RSTARR 

WRE=WR/ROR+UR 

UH(1)=UL-CL 

UW(2)=USTAR-CSTARL 

UW(3)=USTAR 

UW(4)=WRE 

UW(5)=WRE 

GO TO 5 

C THE CASE SS 

2 ITYPE=NCASE 
HELEML=.TRUE. 

HELEMR= .TRUE. 

21 N=N+1 

IF (N.GT.NMAX) GO TO 7003 

ZETAF=ZETA 

PF=ZETAFX5«NG 

PPL=PF/PL - 

PPR=PF/PR 

EVERL=PPL-1 .DO 

EVERR=PPR-1 .DO 

SQRL=DSQRT(1 .D0+G6XEVERL) 

SQRR = DSQRT(1 .D0+G65«EVERR) 
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USL=UL-CLGXEVERL/SQRL 
USR = UR+CRG5(EVERR/SQRR 
DU=USL-USR 

IF (DABS(DU) .LE.UMIDA) GO TO 20 

DUDZL=-NGXCLGX(PPL/ZETAF)X(1 .D0+G9XEVERD/SQRLXX3 
DUDZR= NG3«CRG3<(PPR/2ETAF)X(1 .D0+G95<EVERR)/SQRRXX3 
ZETA=ZETAF+DU/ ( DUDZR-DUDZL ) 

GO TO 21 

0 CONTINUE 
USTAR=(USL+USR)/2.D0 
IFCDABS(USTAR).LT.EPSXUMAX) USTAR=0. 
PSTAR=(PPLXPL+PPRXPR)/2.D0 
WWR=G11X( USTAR-URJXROR 
WR=WWR+DSQRT(GRXX2+NWRX5(2) 

WNL=-G11X(USTAR-UL)XR0L 

WL=WNL+DSQRT(GLXX2+WWLXX2) 

RSTARL=R0LXNL/(WL + R0LX(USTAR-UD) 
RSTARR=R0RXNR/(WR-R0RX(USTAR-UR)) 

GSTARL = DSQRT(GAMA3(PSTAR*RSTARL) 

GSTARR=DSQRT(GAMAXPSTARXRSTARR) 

CSTARL=GSTARL/RSTARL 

CSTARR=GSTARR/RSTARR 

WLE=-WL/ROL+UL 

HRE=WR/ROR+UR 
UN(1)=NLE 
UW(2)=WLE 
UN(3)=USTAR 
UN(9)=WRE 
UW(5)=WRE 
GO TO 5 
THE CASE SE 
ITYPE=NCASE 
HELEML=.TRUE. 

HELEMR=. FALSE. 

1 N=N+1 

IF (N.GT.NMAX) GO TO 7003 
ZETAF=ZETA 

UER = UR+G73<CR3<(ZETAF-ZETAR)/ZETAR 
PPL = (ZETAF/ZETAU3«5(NG 
EVERL=PPL-1 .DO 
SQRL = DSQRT(1 .D0 + G63fEVERL) 

USL=UL-CLGXEVERL/SQRL 

DU=USL-UER 

IF (DABS(DU) .LE.UMIDA) GO TO 30 

DUDZL = -NGXCLGX(PPL/ZETAF)X(1 .D0+G9XEVERL)/SQRLX5<3 
ZETA=ZETAF+DU/ ( DUDZR-DUDZL ) 

GO TO 31 
0 CONTINUE 

USTAR=(USL+UER)/2.D0 
IFCDABS(USTAR) .LT.EPSXUMAX) USTAR=0. 

PSTAR=PPLXPL 

CSTARR=CR-(UR-USTAR)/G7 

RSTARR=GAMAXPSTAR/CSTARRXX2 

GSTARR=CSTARR5(RSTARR 

NWL = -G113<(USTAR-UL)3<ROL 

NL=WWL+DSQRT(GLXX2+HHLXX2) 

WLE=-WL/ROL+UL 

RSTARL=R0LXHL/(HL + R0LX(USTAR-UD) 

GSTARL = DSQRT(GAMA5(PSTARXRSTARL ) 
CSTARL=GSTARL/RSTARL 
UN(1)=WLE 
UH(2)=WLE 
UW(3)=USTAR 
UW(4)=USTAR+CSTARR 
UW(5)=UR+CR 
GO TO 5 
THE CASE EE 
ITYPE=NCASE 
HELEML=. FALSE. 

HELEMR=. FALSE. 

PSTAR=ZETAXXNG 

USTAR=UL-G7XCLX(ZETA-ZETAL)/ZETAL 



CHA19A1 

CHA19A2 

CHA1AA3 

CHA19AA 

CHA1AA5 

CHA19A6 

CHA1A97 

CHA19A8 

CHA1999 

CHA1A50 

CHA1A51 

CHA1A52 

CHA1A53 

CHAIASA 

CHA1A55 

CHA1956 

CHA1A57 

CHA1458 

CHA1A59 

CHA1460 

CHA1A61 

CHA1462 

CHA1463 

CHA1464 

CHA1465 

CHA1466 

CHA1467 

CHA1468 

CHA1469 

CHA1470 

CHA1471 

CHA1472 

CHA1473 

CHA1474 

CHA1475 

CHA1476 

CHA1477 

CHA1478 

CHA1479 

CHA1480 

CHA1481 

CHA1482 

CHA1483 

CHA1484 

CHA1485 

CHA1486 

CHA1487 

CHA1483 

CHA1489 

CHA1490 

CHA1491 

CHA1492 

CHA1493 

CHA1494 

CHA1495 

CHA1496 

CHA1497 

CHA1498 

CHA1499 

CHA1500 

CHA1501 

CHA1502 

CHA1503 

CHA1504 

CHA1505 

CHA1506 

CHA1507 

CHA1508 

CHA1509 

CHA1510 

CHA1511 

CHA1512 



45 



FILE*. 



CHARGER FORTRAN A1 



61 



667 



666 



IF(DABS(USTAR) .LT.EPSxUMAX) USTAR=0. 

CSTARL=CL-KUL-USTAR)/G7 

CSTARR=CR-(UR-USTAR)/G7 

RSTARL=GAMAXPSTAR/CSTARL3(5«2 

RSTARR=GAMA5<PSTAR/CSTARR5«5€2 

GSTARL=RSTARL5«CSTARL 

GSTARR=RSTARRXCSTARR 

UW(1)=UL-CL 

UW(2)=USTAR-CSTARL 

UW(3)=USTAR 

UW(A)=USTAR+CSTARR 

UW(5)=UR+CR 

N=1 

GO TO 5 
CONTINUE 
DO 6 K=l,6 
NFLUX=K 

IF (UW(K) .GE.O.) GO TO 61 

CONTINUE 

NFLUX=6 

CONTINUE 

NC14(NCASE)=NC14(NCASE)*H 

CASEAV(NCASE)=CASEAV(NCASE)-t-DFLOAT(N) 

NF16 INFLUX) =NF16(NFLUX)*H 
IF(NTRY.GE.2)G0 TO 666 
IF(I.NE.2.AND.I.NE.L) GO TO 666 

PRINT 667, I,NFLUX,NCASE,PL,UL,ROL,PR,UR,ROR,USTAR,PSTAR,RSTARL, 

1 RSTARR,(KK,UW(KK),KK=1,6) 

FORMATI/IX, »I,NFLUX,NCASE=',3I5/1X, 'PL,UL,R0L,PR,UR,R0R=»,6D12.4/ 

1 IX, 'USTAR,PSTAR,RSTARL,RSTARR=',4D13.4/ 

2 IX, ’KK,UW(KK)=»,6(I4,2X,D13.4)/) 

NTRY=NTRY+1 



CONTINUE 
RETURN 

7001 CONTINUE 

PRINT 7101, PL,UL,PR,UR,ZETAL,ZETAR,SLL,SRR,NL,NR,I 

7101 F0RMAT(//1X, 'FROM RIEMAN. AN IMPOSSIBLE CASE OF EXPANSION/SHOCK' 

1 //IX, ’PL,UL,PR,UR=',4D25.14// 

2 IX, *ZETAL,ZETAR,SLL,SRR=',4D25.14// 

3 IX, 'NL,NR,I=',3I10//) 

CALL S0FC7001 ') 

7002 CONTINUE 

PRINT 7102, ZETA,DUDZL,DUDZR,ZETAL,ZETAR,PL,UL,PR,UR,N,NCASE,I 

7102 F0RMAT(//1X, ’FROM RIEMAN. NEGATIVE PRESSURE AT THE INTERSECTION' 

1 IX, 'OF L AND R EXPANSION BRANCHES'// 

2 IX, 'IT MEANS THAT A CAVITATION TENDS TO FORM. THIS', 

3 IX, 'POSSIBILITY IS EXCLUDED IN PRESENT VERSION'// 

4 IX, ’ZETA,DUDZL,DUDZR,ZETAL,ZETAR,PL,UL,PR,UR=',9D10.3// 

5 IX, ’N,NCASE,I=',3I10//) 

CALL S0FC7002') 

7003 CONTINUE 

PRINT 7103, I,N,NCASE,DU,UMIDA,EPS,PL,UL,PR,UR, 

1 ZETA,ZETAF,ZETAL,ZETAR,DUDZL,DUDZR 

7103 F0RMATC//1X, 'FROM RIEMAN. NUMBER OF ITERATIONS EXCEEDED.'// 

1 IX, 'I,N,NCASE,DU,UMIDA,EPS=',3I6,3D18.6// 

2 IX, ’PL,UL,PR,UR,ZETA,ZETAF=’,6D18.10// 

3 IX, 'ZETAL,ZETAR,DUDZL,DUDZR=’,4D18.10//) 

CALL SOF( '7003') 

RETURN 

END 

C$OPTIQNS LIST 
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SUBROUTINE MAGACL, I, MIN) h 

IMPLICIT REAL^8(A-H,0-Z,$) 

DIMENSION MIN(L) 

COMMON /GAM/GAMA, NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11 

1 ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23 

2 ,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,G35 
REALMS NG,MU2 

COMMON/DETO/QDET, PCJDET, RCJDET, UCJDET, DCJDET, PODET, ROODET, 
1 RATE,TEMPC 

COMMON /STEPO/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR, 

1 RSTARL,RSTARR,GSTARL,GSTARR, 
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2 CL,CR,CSTARL,CSTARR,SL,SR,WL,WR,UH(6) CHA1588 

3 ,LAMDAL,LAMDAR,RATEL, RATER, TEMPL,TEMPR,TEMPSL,TEMPSR CHA1589 

4 ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR CHA1590 

REALX8 LAMDAL,LAMDAR CHA1591 

LOGICAL HELEML,HELEMR CHA1592 

COMMON /STEP1/DUIDT,DPIDT,DGIDTL,DGIDTR,DRIDTL,DRIDTR CHA1593 

2 ,ASTARL,ASTARR,LAMDSL,LAMDSR,DSDAL,DSDAR,DZDAL,DZDAR CHA1594 

3 ,RAT,SH CHA1595 

4 ,BETACL,BETACR,DSDASL,DSDASR,DZDASL,DZDASR CHA1596 

REALX8 LAMDSL,LAMDSR,DSDAL,DSDAR,DZDAL,DZDAR CHA1597 

COMMON /GRADS/DUDXI L , DPDXI L , DGDXI L , DRDXI L , DZDXI L , DSDXI L , CHA1598 

1 DUDXIR,DPDXIR,DGDXIR,DRDXIR,DZDXIR,DSDXIR CHA1599 

COMMON /AB/A(50) CHA1600 

REAL3<8 LU,LP,LRO,LLAMDA CHA1601 

DATA EPS/1. D-6/ CHA1602 






WE HERE SOLVE FOR THE TIME-DERIVATIVES ALONG THE CONTACT SURFACE 
NAMELY DUIDT,DPIDT. FROM THESE WE ALSO OBTAIN THE OTHER 
TIME-DERIVATIVES (SEE COMMON /STEPl/) . 

WE COMPUTE THE COEFFICIENTS FOR TWO EQUATIONS FOR DUIDT,DPIDT. 
ARE AALXDUIDT+BBLXDPIDT=DDL 

AARXDUIDT+BBRXDPIDT=DDR 



, CHA1604 

CHA1605 
CHA1606 
THESECHA1607 
CHA1608 
CHA1609 






IF(SH.LE.EPS)RAT=0. 

LEFT SIDE OF CONTACT 

IF ( .NOT.HELEML) GO TO 12 

11 CONTINUE 
LEFT SHOCK 

DP=PSTAR-PL 

DU=USTAR-UL 

Z2=0 . 5D0/(PSTAR+MU25(PL ) 

LU = DUX(0.5D0XROL+MU2XZ25(GL5<5<2)-GLXX2/WL-WL 

LRO=-0.5DOXDP/ROL 

LP=-2.D0-MU2XZ2XDP 

AAL=2.D0-Z2XDP 

BBL=Z23(DU+WL/GSTARL3<5(2+1 . DO/WL 
DDL=LUxDUDXIL+LROXDRDXIL+LPXDPDXIL 
DDL=DDL-WL5(USTARXRAT/RSTARL 
1 +ULXRATX(-GAMA3<PL/WL + DUX(GAMAXPLXMU2XZ2+0.5D0)) 

GO TO 10 

12 CONTINUE 
LEFT RAREFACTION 

A1=DUDXIL+DPDXIL/GL 

BETA=GSTARL/GL 

SQB=DSQRT(BETA) 

ASTARL=A1-(CL/(G15XSL))XDSDXILX(BETAXXG5-1 . DO) 

AAL=1 .DO 

BBL=1 .DO/GSTARL 

DDL=-GSTARLXASTARL/SQB 

DSDAL=DSDXIL 

DZDAL=DZDXIL 

DSDASL=DSDXILXSQB 

DZDASL=DZDXILXSQB 

GEOM=RAT^( (GAMA-1 . DO )XUL+2 . DO^CL ) X 
1 (BETAXXG13-1 . DO )/( ROLX( GAMA-3 . DO ) ) 

1 -4. D0XRAT5(CLX(BETAXXG14-1 . DO )/ ( ROLX( 3 . DOXGAMA-5 . DO) ) 

ASTARL=ASTARL-GEOM 
EVER1= GSTARLXGEOM/SQB 
EVER2=-RATXUSTARXCSTARL 
DDL=DDL+EVER1+EVER2 
GO TO 10 
10 CONTINUE 

RIGHT SIDE OF CONTACT 

IF ( .NOT.HELEMR) GO TO 22 
21 CONTINUE 
RIGHT SHOCK 
DP=PSTAR-PR 
DU=USTAR-UR 
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Z2=0.5D0/(PSTAR+MU2XPR) 

LU = DUX(0.5D05«ROR+MU2XZ2XGRX3(2)+GRXx2/WR+WR 

LRO = -0.5D05(DP/ROR 

LP=-2.D0-MU2XZ2XDP 

AAR=2.D0-Z2XDP 

BBR=Z2XDU-WR/GSTARRXX2-1 .DO/WR 
DDR = LU5€DUDXIR+LROxDRDXIR+LP>(DPDXIR 
DDR = DDR+WR5(USTAR3(RAT/RSTARR 
1 +URXRATX(GAMAxpR/WR+DU5f(GAMAxPR3(MU2XZ2+0 .500) ) 

GO TO 20 
22 CONTINUE 
C RIGHT RAREFACTION 

A1=DUDXIR-DPDXIR/GR 

BETA=GSTARR/GR 

SQB=DSQRT(BETA) 

ASTARR = A1 + (CR/(G153«SR))XDSDXIRX(BETAXXG5-1.D0) 

AAR=1.D0 

BBR=-1 .DO/GSTARR 

DDR=GSTARRXASTARR/SQB 

DSDAR=DSDXIR 

DZDAR=DZDXIR 

DSDASR=DSDXIRX$QB 

DZDASR=DZDXIRXSQB 

GE0M=RATX(-(GAMA-1 . D0)*UR+2 . D0XCR)5f(BETA*XG13-l . DO) 

1 /(RORX(GAMA-3.DO)) 

2 -4.DO>6RATXCR5((BETAXxG1A-1.DO)/(RORX(3.DOxGAMA-5.DO)) 
ASTARR=ASTARR+GEOM 

EVER1=GSTARRXGE0M/SQB 
EVER2=RATXUSTARXCSTARR 
DDR=DDR+EVER1+EVER2 
GO TO 20 
20 CONTINUE 

DET=AAL>(BBR-AARXBBL 
DUIDT = (DDL5tBBR-DDR3(BBL)/DET 
DPIDT=-(DDLXAAR-DDRXAAL)/DET 
DRIDTL = DPIDT/CSTARLX5(2 
DRI DTR = DPIDT/CSTARR5(X2 
RETURN 
END 



SUBROUTINE FLUXE( L , I , MIN) 

IMPLICIT REAL5<8(A-H,0-Z,$) 

DIMENSION MIN(L) 

COMMON /AB/A(50) 

EQUIVALENCE ( DT, A( 4) ) , ( NCYC, A( 12) ) 

COMMON /GAM/ GAMA , NG, MUZ , G1 , G2 , G3 , G4 , G5 , G6 , G7 , G8 , G9 , G1 0 , G1 1 

1 ,G12,G13,G14,G15,G16,G17,G18,G19,G20,G21,G22,G23 

2 ,G24,G25,G26,G27,G28,G29,G30,G31,G32,G33,G34,G35 
REAL)C8 NG,MU2 

COMMON /GRADS/DUDXI L , DPDXI L , DGDXI L , DRDXI L , DZDXI L , DSDXI L , 

1 DUDXIR, DPDXIR, DGDXI R,DRDXIR,DZDXIR, DSDXI R 

COMMON /STEPO/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR, 
RSTARL,RSTARR,GSTARL,GSTARR, 

CL , CR, CSTARL , CSTARR, SL , SR, ML , MR, UW( 6 ) 
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CHA1700 
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CHA1703 
CHA1704 
CHA1705 
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CHA1711 
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,LAMDAL,LAMDAR,RATEL, RATER, TEMPL,TEMPR,TEMPSL,TEMPSR CHA1713 






1 

2 

3 

4 



, ZL , ZR, ZSTARL , ZSTARR, NFLUX , HEL EML , HEL EMR CHA17 1 4 

REALMS LAMDAL, LAMDAR CHA1715 

LOGICAL HELEML,HELEMR CHA1716 

COMMON /STEP1/DUIDT,DPIDT,DGIDTL,DGIDTR,DRIDTL,DRIDTR CHA1717 

2 ,ASTARL,ASTARR,LAMDSL,LAMDSR,DSDAL,DSDAR,DZDAL,DZDAR CHA1718 

3 ,RAT,SH CHA1719 

4 ,BETACL,BETACR,DSDASL,DSDASR,DZDASL,DZDASR CHA1720 

REAL5«8 LAMDSL,LAMDSR,DSDAL,DSDAR,DZDAL,DZDAR CHA1721 

COMMON/DETO/QDET,PCJDET,RCJDET,UCJDET,DCJDET,PODET,ROODET, CHA1722 

1 RATE,TEMPC CHA1723 

COMMON /FI/FIH1,FIH2,FIH3,UXN,PXN,GXN,R0XN,ZXN CHA1724 

1 ,GIH CHA1725 

2 ,FIH4,ZMD0TL,ZMD0TR CHA1726 

REAL5(8 LAMDAO CHA1727 

C RO,U,P,Z AND THEIR (XI, T) DERIVATIVES AT EULERIAN POINT X=X(I). CHA1729 

DT2=DT/2.D0 CHA1731 
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GO TO (1,2,3,A,5,6),NFLUX 
CONTINUE 

NFLUX=1. LINE X=0 IS TO THE LEFT OF LEFT WAVE. 

UX=UL 
PX=PL 
ROX=ROL 
ZX = ZL 
GX=GL 

DUDXIX=DUDXIL 

DPDXIX=DPDXIL 

DRDXIX=DRDXIL 

DZDXIX=DZDXIL 

DUDTX=-DPDXIL 

DR0DTX=-R0LX5<25<DUDXIL 

DPDTX=-GLX3<2xdUDXIL 

DRODTX=DRODTX-RATXROLXUL 

DPDTX=DR0DTX3«CLX3<2 

DZDTX=0. 

GO TO 9 
CONTINUE 

NFLUX=6. LINE X=0 IS TO THE RIGHT OF RIGHT WAVE. 

UX=UR 

PX=PR 

ROX=ROR 

ZX=ZR 

GX=GR 

DUDXIX=DUDXIR 

DPDXIX=DPDXIR 

DRDXIX=DRDXIR 

DZDXIX=DZDXIR 

DUDTX=-DPDXIR 

DPDTX = -GR5«3(2^DUDXIR 

DRODTX = -ROR5«5(23<DUDXIR 

DRODTX=DRODTX-RAT3«ROR5«UR 

DPDTX=DRODTX5€CR5«5€2 

DZDTX=0. 

GO TO 9 
CONTINUE 

NFLUX=2. SONIC CASE (LEFT). 

BETA0 = (MU25«(UL/CL+G7) .D0/MU2) 

SQBO=DSQRT(BETAO) 

A1=DUDXIL+DPDXIL/GL 

A0 = A1-(CL/(G155<SL))5(DSDXIL3«(BETA05«^G5-1 . DO) 

EVER1=-((GAMA-1 .D0)5<UL + 2.D05(CL)5((BETA05(3<G13-1 .D0)/(GAMA-3.D0) 
EVER2 = A. D0S«CLX(BETA05«3<G1A-1 .D0)/(3.D05(GAMA-5.D0) 

EVER=( EVER1 + EVER2)5(RAT/R0L 
AO=(AO+EVER) 

DPDAX=GL5<BETA05«A0 
C0 = MU2s<(UL+G75(CL ) 

IFCCO.LT.O.) CALL SOFCFLUXE 2. CO NEGATIVE.') 

UX=C0 

ROX = GL5(BETAO/CO 
ZX = ZL 

PX=ROX5(C05«5(2/GAMA 

GX=ROX5^CO 

DPDAX = DPDAX+RAT5<UX5<C05(SQB0 
DUDBX = -CL5«BETA05<?((-1 .D0/G4)/GA 
DPDBX=PL5(BETA05«5<MU2/G6 
DRODBX = ROL5<BETA05«5<(-MU2)/G‘4 
DSDAX = SQB05(DSDAL 
DZDAX=SQBO^DZDAL 

DRODAX=DPDAX/CO^?«2-(ROX/(GAMA?«SL) )5€DSDAX 
DUDAX=A0 

DGDAX = 0 .5D05eGAMA5((PX5«DR0DAX+R0X5«DPDAX)/GX 

GO TO 9 

CONTINUE 
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NFLUX=5. SONIC CASE (RIGHT). 

BETA0 = (MU2X(-UR/CR+G7))X5(<1.D0/MU2) 

SQBO=DSQRT(BETAO) 

A1=DUDXIR-DPDXIR/GR 

A0=A1 + (CR/(G153«SR))xDSDXIR5((BETA03(XG5-1 .DO) 

EVERl =( -(GAMA-1 .DO) XUR+2.D0XCR)X(BETA0XXG13-1 .DO)/ (GAMA-3. DO) 
EVER2 = -A.D0XCRX(BETA0X5€G1A-1 .D0)/(3.D0XGAMA-5.D0) 
EVER=(EVER1+EVER2)XRAT/R0R 
A0=(A0+EVER) 

DPDAX=-GR3(BETA03(A0 

C0=MU2X(-UR+G7XCR) 

IF(CO.LT.O.) CALL SOFCFLUXE 5. CO NEGATIVE.') 

UX=-C0 

R0X=GRXBETA0/C0 

ZX=ZR 

PX=ROX3(COXX2/GAMA 

Gx=R0X5^C0 

DPDAX=DPDAX-RATXUXXC0XDSQRT(BETA0) 

DUDBX=CRXBETA0XX(-1 . D0/GA)/GA 

DPDBX=PRXBETA0XXMU2/G6 

DR0DBX=R0RXBETA0X5((-MU2)/GA 

DSDAX=SQBOXDSDAR 

DZDAX=SQB03«DZDAR 

DRODAX=DPDAX/COXX2-(ROX/(GAMAXSR))XDSDAX 

DUDAX=A0 

DGDAX=0 . 5DOXGAMAX(PXXDRODAX+ROXXDPDAX)/GX 
GO TO 9 

3 CONTINUE 

NFLUX=3. LINE X=0 IS BETWEEN THE LEFT WAVE AND THE CONTACT. 

UX=USTAR 

PX=PSTAR 

ROX=RSTARL 

ZX=ZL 

GX=GSTARL 

DUDXIX=-DPIDT/GSTARL 3(5(2 
DPDXIX=-DUIDT 

DUDXIX = DUDXIX-RAT3«USTAR/RSTARL 

DZDXIX=DZDXIL 

DZDTX=0. 

IF ( .NOT.HELEML) GO TO 32 

31 CONTINUE 
C LEFT SHOCK. 

DRDXIX=(RSTARL/WL)3(3(23((3.D03(DUIDT 

1 +DPIDT3((1 .D0+3.D03((WL/GSTARL)3(3(2)/WL 

2 +DUDXIL5(WL3(((GL/WL)3(X2+3.D0)+3.D03(DPDXIL 

3 +DRDXIL3((WL/R0L)3(3(2) 

EV ERl =UL 3(RSTARL3( 3(2XRAT3(( (GL/WL)3( 3(2+1 .DO)/ (ROL3(WL) 

EVER2=2 . D03(RSTARL3(USTAR3(RAT/WL 
DRDXIX=DRDXIX+EVER1+EVER2 
DRODTX=-DUDXIX3(ROX3(3(2 
GO TO 33 

32 CONTINUE 
BETA=GSTARL/GL 
SQB=DSQRT(BETA) 

DPDA=ASTARL3(GSTARL 

DPDA=GSTARL3((ASTARL + RAT3(USTAR3(CSTARL/(GL3( SQB)) 

GA1=1 .D0/GA+0.5D0 

DRODA = (DRDXIL-DPDXIL/(CL3(CL)) 3(BETA3(3(GA1 + DPDA/ (CSTARL3(3(2) 

DRDXIX= DR0DA/SQB+DPIDT/(GSTARL3(CSTARL3(3(2) 

DRODA = DPDA/CSTARL3(;(2-(RSTARL/(GAMA3(SL))3(DSDASL 

DR0DTX=-DUDXIX3(R0X3(3(2 

DRDXIX=DRODA/SQB+DRODTX/GSTARL 

33 CONTINUE 
DUDTX=DUIDT 
DPDTX=DPIDT 
GO TO 9 

4 CONTINUE 
C 
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CHA1828 

CHA1829 

CHA1830 

CHA1831 

CHA1832 
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CHA1837 

CHA1838 

CHA1839 

CHA1840 

CHA1841 

CHA1842 

CHA1843 

CHA1844 

CHA1845 

CHA1846 

CHA1847 

CHA1848 

CHA1849 

CHA1850 

CHA1851 

CHA1852 

CHA1853 

CHA1854 

CHA1855 

CHA1856 

CHA1857 

CHA1858 

CHA1859 

CHA1860 

CHA1861 

CHA1862 

CHA1863 
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CHA1874 

CHA1875 
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GEP FORTRAN A1 



PAG 



NFLUX=4. LINE X=0 IS BETWEEN THE CONTACT AND THE RIGHT WAVE. 

DPDXIX=-DUIDT 

UX=USTAR 

PX=PSTAR 

ROX=RSTARR 

ZX=ZR 

GX=GSTARR 

DUDXIX=-DPIDT/GSTARR5(S<2 

DUDXIX=DUDXIX-RAT5«USTAR/RSTARR 

DPDXIX=-DUIDT 

DZDXIX=DZDXIL 

DZDTX=0 . 

IF ( .NOT.HELEMR) GO TO 42 

41 CONTINUE 
RIGHT SHOCK 

DRDXIX=(RSTARR/WR)5(3(23((3.3(DUIDT 

1 -DPIDT5<(1 .D0+3.D0X(WR/GSTARR)S<X2)/WR 

2 -DUDXIR3«WRX( (GR/WR)X3(2+3. D0) + 3.D05<DPDXIR 

3 +DRDXIRX(WR/R0R)5<9(2) 

EVER1=URxRSTARRxx2J(RAT3(( (GR/WR)5(J(2+1 .DO)/(ROR5(WR) 

EVER2=2 . D0XRSTARR5<USTAR5(RAT/WR 

DRDXIX=DRDXIX-EVER1-EVER2 

DRODTX=-DUDXIX5(ROX5(3(2 

GO TO 43 

42 CONTINUE 
RIGHT RAREFACTION 

BETA=GSTARR/GR 

SQB=DSQRT(BETA) 

DPDA=-ASTARRXGSTARR 

DPDA=-GSTARR3(( ASTARR+RAT5<USTAR5eCSTARR/(GR5( SQB ) ) 
G41=1.D0/G4+0.5D0 

DRODA=(DRDXIR-DPDXIR/(CR3«CR)) 3(BETA5(5(G41 + DPDA/ (CSTARR5(3<2 ) 

DRDXIX= DR0DA/SQB-DPIDT/(GSTARR5<CSTARRX)(2) 

DRODA = DPDA/CSTARR5(5(2-(RSTARR/(GAMA3(SR))5«DSDASR 
DRODTX = -DUDXIX3(ROX5(5(2 
DRDXIX=DRODA/SQB-DRODTX/GSTARR 

43 CONTINUE 
DUDTX=DUIDT 
DPDTX=DPIDT 
GO TO 9 

9 CONTINUE 
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CHA1913 

CHA1914 

CHA1915 
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FLUXES CENTERED AT TIME T(N+l/2) AT EULERIAN POINT X=X( I) . CHA1919 






FIl=ROX5(UX 

FI2 = ROXMUX5(K2+PX 

FI2=FI2-PX 

FI3 = UX5((G12?<PX+0 .5D05<ROX>(UXX5(2) 

FI4 = ZX5«ROX5(UX 
FI3 = FI3 + QDET5(FI4 
ROUOO = ROX)(UX 

GO T0(10,20,30,40,50,60), NFLUX 
10 CONTINUE 
60 CONTINUE 

DFDXI1 = DRDXIX5(UX+R0X3(DUDXIX 

DFDXI2 = DRDXIX5(UX5(^(2 + 2.DOXROX5(UX5(DUDXIX+DPDXIX 
DFDXI2=DFDXI2-DPDXIX 

DFDXI3 = DUDXIXK(G125(PX+0.5D05(ROX5(UX3(3(2) 

1 +UX5((G125(DPDXIX+0.5D03(DRDXIX5(UX5(3(2+ROXxUX5(DUDXIX) 

DFDXI4 = ZX5<DFDXI1 + R0X5(UX5(DZDXIX 
DFDXI3=DFDXI3+QDETmDFDXI 4 
DFIDTl = DRODTX3(UX+ROX3<DUDTX 

DFIDT2=DRODTX5(UX5(K2+2.D05(ROX5(UX5(DUDTX+DPDTX 

DFIDT2=DFIDT2-DPDTX 

DFIDT3 = DUDTX5((G12¥PX+0 .5D03«ROX3(UX5(5(2) 

1 +UX3((G125(DPDTX + 0.5D03(DRODTX5(UX3(5(2+ROX5(UX3(DUDTX) 

DFIDT4 = ZX5(DFIDT1 + R0X^ZXMDZDTX 
DFIDT3 = DFIDT3 + QDET5«DFIDT4 
FID0T1=-R0U00xdFDXI1+DFIDT1 
FIDOT2 = -ROU003(DFDXI2 + DFIDT2 
FIDOT3 = -ROU005«DFDXI3 + DFIDT3 
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FILE* CHARGER FORTRAN A1 





FIDOT4=-ROUOOXDFDXI4+DFIDTA 


CHA1948 




UXDOT=-ROUOOXDUDXIX+DUDTX 


CHA1949 




PXDOT=-ROUOOXDPDXIX+DPDTX 


CHA1950 




ROXDOT=-ROU005^DRDXIX+DRODTX 


CHA1951 




2XDOT = -ROU003(DZDXIX+DZDTX 


CHA1952 




FIH1=FI1+DT2XFID0T1 


CHA1953 




FIH2=FI2+DT2XFID0T2 


CHA1954 




GIH=PX+DT2XPXD0T 


CHA1955 




FIH3=FI3+DT2XFID0T3 


CHA1956 




FIH4=FIA+DT2XFID0T4 


CHA1957 




UXN=UX+DTXUXDOT 


CHA1958 




PXN=PX+DTXPXD0T 


CHA1959 




ROXN=ROX+DTXROXDOT 


CHA1960 




ZXN=ZX+DTXZXD0T 


CHA1961 




IF(ZXN.LT.O. ) ZXN=0. 


CHA1962 




GO TO 90 


CHA1963 


20 


CONTINUE 


CHA1964 




EV0=GL5(DSQRT(BETA0) 


CHA1965 


201 


CONTINUE 


CHA1966 




DFIDA1=DRODAX5(UX+ROX3«DUDAX 


CHA1967 




DFIDA2=DRODAXXUX5(X2+2 . DOXROX^UXXDUDAX+DPDAX 


CHA1968 




DFIDA2=DFIDA2-DPDAX 


CHA1969 




DFIDA3 = DUDAXX(G12XPX+0 .5D05(ROX5(UX5«X2) 


CHA1970 




1 +UXX(G12XDPDAX+0 . 5DOXDRODAX3<UXXX2+ROXXUXXDUDAX) 


CHA1971 




DFIDA4=ZX3<DFIDA1+R0X>(UX5(DZDAX 


CHA1972 




FIDOT1=-EVOXDFIDA1 


CHA1973 




FIDOT2=-EVOXDFIDA2 


CHA1974 




FIDOT3=-EVOXDFIDA3 


CHA1975 




FIDOT4=-EV05(DFIDA4 


CHA1976 




FIH1=FI1+DT2XFID0T1 


CHA1977 




FIH2=FI2+DT2XFID0T2 


CHA1978 




FIH3=FI3+DT2XFID0T3 


CHA1979 




FIH4=FI4+DT2XFID0T4 


CHA1980 




GA=DGDAX 


CHA1981 




IF(NFLUX.EQ.5)GA=-GA 


CHA1982 




droua=uxxdrodax+roxxdudax 


CHA1983 




BETAPR=0.5D036DSQRT(BETA0)x(GA-DROUA) 


CHA1984 




FIH2=FIH2-DPDBXXBETAPRXDT2 


CHA1985 




UXDOT = -EVO^DUDAX+BETAPR5«DUDBX 


CHA1986 




PXDOT=-EV03<DPDAX+BETAPR5<DPDBX 


CHA1987 




GIH = PX+DT25(PXDOT 


CHA1988 




roxdot=-evoxdrodax+betaprxdrodbx 


CHA1989 




ZXD0T = -EV05«DZDAX 


CHA1990 




UXN=UX+DTXUXDOT 


CHA1991 




PXN=PX+DTXPXD0T 


CHA1992 




ROXN=ROX+DTXROXDOT 


CHA1993 




ZXN=ZX+DTXROXDOT 


CHA1994 




IFCZXN.lt. 0.) ZXN=0. 


CHA1995 




GO TO 90 


CHA1996 


50 


CONTINUE 


CHA1997 




EVO=-GRXDSQRT(BETAO) 


CHA1998 




GO TO 201 


CHA1999 


30 


CONTINUE 


CHA2000 


40 


CONTINUE 


CHA2001 




GO TO 60 


CHA2002 


90 


CONTINUE 


CHA2003 




RETURN 


CHA2004 




END 


CHA2005 
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Figure A- 1. Piecewise Linear Distribution of Flow Variables in Cells 
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Figure A-2. Intersection of Right and Left Adiabats for Solving Riemann Problem 
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Figure A-3. Wave Diagram Representing Solution to Riemann Problem 
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APPENDIX B. Code for Re-Normalizing the Air Impulse 



1 IMPLICIT REALX8(A-H,0-Z) 

C CODE RENORM — C TRANSFORMATION OF TOTAL REFLECTED IMPULSE FROM 
C BAKER'S CHART TO SPACE-NORMALIZED VALUES. 

C DATA FROM FIG. 6.3 (SUPPLEMENT) IN BAKER'S BOOK "EXPLOSIONS IN AIR" 

2 REALXA RB,IB,RS,IS,ISBARE 

3 DIMENSION RB( 21 ) , IB( 21 ) 

4 DIMENSION RS( 21 ) , IS(21 ) , ISBARE(21 ) 

5 DATA RB/. 05, .06, .07, .08, .09, .1, .2, .3, .4, .5, .6, .7, .8, .9,1., 

1 2. ,3. ,4. ,5. ,6. ,7./ 

6 DATA IB/4.4,3.06,2.30,1.83,1.50,1.27, .457, . 293, . 221 , . 178 , . 149 , 

1 .128, .113, .099, .0885, .0376, .0236, . 0173, . 0136, . 0113, .0095/ 

7 PAI=4.D0*DATAN(1.D0) 

8 G=1.4D0 

9 PA=0.1D0 

10 RHOA=1.3DO 

11 RHO0=1800.D0 

12 Q0=4.D0 

13 BETA=DSQRT(RHOA/RHOO)*(PA/(RHOOXQO))xx(1.DO/6.DO) 

14 GOREM=(3.DO/DSQRT(2.DOXG))X(4.DOXPAI/3.DO)XX(1.DO/3.DO) 

15 BETA=BETAxGOREM 

16 DELTA=( (4.DOXPAI/3.DO)X(RHOOXQO/PA) ) XX( 1 . DO/3 . DO ) 

17 PRINT 11, BETA, DELTA 

18 11 F0RMAT(/1X, 'RESULTS WITH BETA, DELTA= ' , 2D16 . 7// 

1 IX,' N',' RB ',' IB ',2X, 

2 ' RS ',' IS ',2X,' ISBARE '/) 

19 DO 1 N=l,21 

20 RS(N)=RB(N)XDELTA 

21 IS(N)=IB(N)XBETA 

22 ISBARE(N)=1 .D0/RS(N)XX2 

23 PRINT 2, N,RB(N),IB(N),RS(N),IS(N),ISBARE(N) 

24 2 F0RMAT(1X,I4,2E12.4,2X,2E12.4,2X,E12.4) 

25 1 CONTINUE 

26 END 



RESULTS WITH BETA 


,DELTA= 0. 


.1204163D-01 


0.6706157D+02 




N 


RB 


IB 


RS 


IS 


ISBARE 


1 


0.5000E-01 


0. AAOOE+01 


0.3353E+01 


0.5298E-01 


0.8894E-01 


2 


0.6000E-01 


0.3060E+01 


0 .4024E+01 


0.3685E-01 


0.6177E-01 


3 


0.7000E-01 


0.2300E+01 


0.4694E+01 


0.2770E-01 


0.4538E-01 


4 


0.8000E-01 


0.1830E+01 


0.5365E+01 


0.2204E-01 


0.3474E-01 


5 


0.9000E-01 


0.1500E+01 


0 .6036E+01 


0.1806E-01 


0.2745E-01 


6 


O.lOOOE+00 


0.1270E+01 


0.6706E+01 


0.1529E-01 


0.2224E-01 


7 


0.2000E+00 


0. A570E+00 


0.1341E+02 


0.5503E-02 


0.5559E-02 


8 


0 .3000E+00 


0.2930E+00 


0.2012E+02 


0.3528E-02 


0.2471E-02 


9 


0 .4000E+00 


0.2210E+00 


0.2682E+02 


0.2661E-02 


0.1390E-02 


10 


0.5000E+00 


0.1780E+00 


0.3353E+02 


0.2143E-02 


0 .8894E-03 


11 


0.6000E+00 


0.1A90E+00 


0.4024E+02 


0.1794E-02 


0.6177E-03 


12 


0.7000E+00 


0.1280E+00 


0.4694E+02 


0.1541E-02 


0 .4538E-03 


13 


0.8000E+00 


0.1130E+00 


0 .5365E+02 


0.1361E-02 


0.3474E-03 


14 


0.9000E+00 


0.9900E-01 


0.6036E+02 


0.1192E-02 


0.2745E-03 


15 


O.lOOOE+01 


0.8850E-01 


0.6706E+02 


0.1066E-02 


0.2224E-03 


16 


0.2000E+01 


0.3760E-01 


0.1341E+03 


0.4528E-03 


0 .5559E-04 


17 


0.3000E+01 


0 .2360E-01 


0.2012E+03 


0 .2842E-03 


0.2471E-04 


18 


0.<4000E+01 


0.1730E-01 


0.2682E+03 


0 .2083E-03 


0.1390E-04 


19 


0.5000E+01 


0. 1360E-01 


0.3353E+03 


0.1638E-03 


0.8894E-05 


20 


0.6000E+01 


0.1130E-01 


0.4024E+03 


0.1361E-03 


0.6177E-05 


21 


0.7000E+01 


0,9500E-02 


0 .4694E+03 


0.1144E-03 


0.4538E-05 
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RENOO 
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RENOO 
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